Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

vtk output for solutions #200

Open
XuliaDS opened this issue Nov 17, 2023 · 14 comments
Open

vtk output for solutions #200

XuliaDS opened this issue Nov 17, 2023 · 14 comments

Comments

@XuliaDS
Copy link

XuliaDS commented Nov 17, 2023

Hello,
I am trying to output the solution in vtk but the file format doesn't seem to match.
I run the following:

mesh.PrintVTK(os.path.join(output_dir, "final_mesh.vtk"), 8)
mesh.Print(os.path.join(output_dir, "final_mesh.mesh"), 8)
x.Save(os.path.join(output_dir, "displacement.gf"))
x.SaveVTK(os.path.join(output_dir, "displacement.vtk"), 'U', 8)

The mesh output is perfect. It has the elements and the materials ID. the solution output (displacement.vtk) does not seem a vtk file. It is just a list of entries.

image

Is this issue a bug or am I doing something wrong ?

Thank you !

@v-dobrev
Copy link
Member

I think SaveVTK is meant to append its output to the same file as PrintVTK -- in VTK the mesh and solution go to the same file. See the documentation here: https://github.com/mfem/mfem/blob/8d75a62e6534be9ca01c7aa15e849dbc7c2a9201/fem/gridfunc.hpp#L749-L752.

@XuliaDS
Copy link
Author

XuliaDS commented Nov 21, 2023

Thanks for the info and it actually makes a lot of sense... (vtk mesh and field altogether).
If I understood correctly, I should do:
ref = 4
mesh.PrintVTK(os.path.join(output_dir, "final_mesh.vtk"), ref)
x.SaveVTK(os.path.join(output_dir, "final_mesh.vtk"), 'U', ref)

I've tried this but unfortunately, when I call SaveVTK from the solution, it overwrites the file rather than appending it.
Any further thoughts ?

@v-dobrev
Copy link
Member

It looks like the issue is in the python wrapper -- in the C++ code the first argument to both PrintVTK and SaveVTK is an output stream which should be the same file -- both methods just append to the stream. In your code snippet, it looks like the python wrapper replaces the first argument with just a file name and internally re-opens the file in SaveVTK and that truncates it and overwrites the output from PrintVTK. Maybe there is a version of SaveVTK in the python version that uses a file handle (or something like that) that can be used in both PrintVTK and SaveVTK to prevent the overwriting? @sshiraiwa, is this possible in the python version?

@sshiraiwa
Copy link
Member

@XuliaDS and @v-dobrev
Understood. The current implementation always creates a new file when the filename is
given as text. I suppose that a work-around is to use StringIO to collect the text data first and save it later.

   from io import StringIO                                                
    output = StringIO()                                                
    output.precesion = 8                                                  
    mesh.PrintVTK(output, 4)                                             
    x.SaveVTK(output, 'U', 4)                                                 
    fid = open("vtkdata", "w")                                              
    fid.write(output.getvalue())                                            
    fid.close()

@v-dobrev I didn't realize that mfem::ostream is used in such a way to append the data for the VTK case. It there any other C++ methods which I should be careful with? Perhaps, in the future, it is better to change the wrapper so that it opens a file in the append mode if the file already exists. But I would rather do this only for the methods that are supposed to append contents, like this one.

@v-dobrev
Copy link
Member

@sshiraiwa, I'm not 100% sure but I think SaveVTK is the only method that is meant to append to a file.

Something similar happens with the socket streams when we send data to GLVis (i.e. the mesh is written to the socket stream and then the grid-function is appended to the same socket stream, typically, using the MFEM-native format with Mesh::Print and GridFunction::Save) but in that case there are no files created.

@XuliaDS
Copy link
Author

XuliaDS commented Dec 1, 2023

Hi!
thank you very much @sshiraiwa and @v-dobrev. I can confirm that this approach works !!
I have added x and the stress values to a VTK and both display perfect. Just for completeness, I write it here again without typo (precesion -> precision)

from io import StringIO
output = StringIO()
output.precision = 8
mesh.PrintVTK(output, 4)
x.SaveVTK(output, 'U', 4)
fid = open("vtkdata", "w")
fid.write(output.getvalue())
fid.close()

many thanks!
Xulia

@XuliaDS XuliaDS closed this as completed Dec 1, 2023
@XuliaDS XuliaDS reopened this Jan 9, 2024
@XuliaDS
Copy link
Author

XuliaDS commented Jan 9, 2024

Hi again,
Extra comment / doubt:
I have tried saving exactly the same mesh but the output vtk from mfem changes the original data.
The image shows: left the original cube (cube.vtk) and right: the converted cube in mfem (converted_mesh.vtk)

image

I run the following:


gmsh_mesh = "cube.msh"
mesh = mfem.Mesh(gmsh_mesh)
output1 = io.StringIO()
output1.precision = 8
refm = 0
mesh.PrintVTK(output1, refm)
fid = open("converted_mesh.vtk", "w")
fid.write(output1.getvalue())


As you can see, it is like every tet has now a surface patch. Is there a way to preserve the orginal mesh ? This extra layers make it very difficult to see inside meshes. For example, if you have two embedded objects and want to look at the interface.

The full script is here:


``
import mfem.ser as mfem
import gmsh
import sys
import os
import io
create_gmsh = True
gmsh_mesh = "cube.msh"
if create_gmsh:

gmsh.initialize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)	
gmsh.model.add("t1")
lc = 0.5
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(1, 0, 0, lc, 2)
gmsh.model.geo.addPoint(1, 1, 0, lc, 3)
gmsh.model.geo.addPoint(0, 1, 0, lc, 4)

gmsh.model.geo.addPoint(0, 0, 1, lc, 5)
gmsh.model.geo.addPoint(1, 0, 1, lc, 6)
gmsh.model.geo.addPoint(1, 1, 1, lc, 7)
gmsh.model.geo.addPoint(0, 1, 1, lc, 8)

gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)

gmsh.model.geo.addLine(5, 6, 5)
gmsh.model.geo.addLine(6, 7, 6)
gmsh.model.geo.addLine(7, 8, 7)
gmsh.model.geo.addLine(8, 5, 8)

gmsh.model.geo.addLine(1, 5, 9)
gmsh.model.geo.addLine(2, 6, 10)

gmsh.model.geo.addLine(3, 7, 11)
gmsh.model.geo.addLine(4, 8, 12)			

gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)  # bottom 
gmsh.model.geo.addCurveLoop([5, 6, 7, 8], 2)  # top 
gmsh.model.geo.addCurveLoop([9, -8, -12, 4], 3)  # left 
gmsh.model.geo.addCurveLoop([1, 10, -5, -9], 4)  # front 
gmsh.model.geo.addCurveLoop([2, 11, -6, -10], 5)  # right
gmsh.model.geo.addCurveLoop([-3, 11, 7, -12], 6)  # back

for j in range(6):
	sj = gmsh.model.geo.addPlaneSurface([j+1], j+1)
loop = gmsh.model.geo.addSurfaceLoop([1,2,3,4,5,6])
vol = gmsh.model.geo.addVolume([loop])
gmsh.model.geo.synchronize()
for j in range(6):
	gmsh.model.addPhysicalGroup(2, [sj], name="box_"+str(j+1))

gmsh.model.addPhysicalGroup(3, [vol], name="volume", tag=1)		

# We can then generate a 2D mesh...
gmsh.model.mesh.generate(3)
gmsh.write(gmsh_mesh)
gmsh.write(gmsh_mesh.replace(".msh", ".vtk"))
gmsh.finalize()

mesh = mfem.Mesh(gmsh_mesh)

output1 = io.StringIO()
output1.precision = 8
refm = 0
mesh.PrintVTK(output1, refm)
fid = open("converted_mesh.vtk", "w")
fid.write(output1.getvalue())
``


@sshiraiwa
Copy link
Member

Comparing mesh.vtk and converted_mesh.vtk from your script. even the number of vertices differs.
The converted one record 400 points, instead of 45, which is 4 x 100, where 100 is the number of elements.
Thus, I suspect that the MFEM VTK output treats all elements separately and does not consider duplication of
vertices. Face information may be treated in a similar manner. I don't know if there is an option to export the vtk in the way you are looking for. @v-dobrev , any suggestion?

@v-dobrev
Copy link
Member

There are two versions of PrintVTK: https://github.com/mfem/mfem/blob/d5aa18f7625032d35e105381f2d34b4a2ffb2db1/mesh/mesh.hpp#L2166-L2174

It looks like you are using the second one -- try using the first one. The second version is for visualization of high-order curved meshes and it breaks down the mesh to individual elements, so that it can support discontinuous, H(div)-, and H(curl)-conforming fields. The first version has the limitation that it only supports linear and quadratic meshes.

@XuliaDS
Copy link
Author

XuliaDS commented Jan 12, 2024

Cool. So the answer was: call the function as simple as possible :)

Ok, this works for the mesh but then how do you match it with the solution ? For example, using example2 in the tutorials (linear elasticity): lets say I have a box and I apply a compression on the top lid and I want to visualize the displacement, then I need to output my solution x.PrintVTK. According to the documentation:
https://github.com/mfem/mfem/blob/2b27d4815f6e549ffb01065e9ad2d3549706ccec/fem/gridfunc.hpp#L755

/** @brief Write the GridFunction in VTK format. Note that Mesh::PrintVTK
must be called first. The parameter ref > 0 must match the one used in
Mesh::PrintVTK.
*/

This means that we need to use the version for mesh VTK with ref parameter. I thought initially, that if I used the refinement = 1, it wouldn'd do any subdivision. I even tried with ref = 0 but it behaves as ref = 1. Then with ref = 2, you get a subdivision and so on.

I tried "fooling" the function by adding a dumb ref:


output = io.StringIO()
output.precision = 8
mesh.PrintVTK(output)

ref = 1
x.SaveVTK(output, 'U', ref)
fid = open("converted_mesh.vtk", "w")
fid.write(output.getvalue())

but then it maps the data to the mesh wrongly. I also get an error on load (shown in screenshot). Which is not surprising since I am giving it a "random" reference

image

I've also noticed that if you give the mesh with reference parameter, the x data saves by point-values in the cells. If you give it without it, you get cell values.

image

An updated version of the script including the MFEM call:

import mfem.ser as mfem
import gmsh
import sys
import os
import io


create_gmsh = True
create_mfem = True
gmsh_mesh = "cube.msh"
lc = 0.1
if create_gmsh:
	gmsh.initialize()
	gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)	
	gmsh.model.add("gmsh to mfem VTK")
	iterX = [0,1,1,0]
	iterY = [0,0,1,1]
	cA = [1, 2, 3, 4, 1]
	pID = 1 ; eID = 1 ; cID = 1; sID = 1; 
	a6 = range(6)
	b_names = {1: "bottom", 2: "top", 3: "left", 4: "front", 5: "right", 6: "back"}
	for zi in range(2):
		for k in range(4):
			print(" ADD POINT ",iterX[k], iterY[k], zi, lc, pID)
			gmsh.model.geo.addPoint(iterX[k],iterY[k],zi, lc, pID)
			pID +=1
	gmsh.model.geo.synchronize()
	for zi in range(2):
		offL = zi * 4
		for k in range(len(cA)-1):
			print(" ADD LINE ", cA[k] + offL, cA[k+1]+ offL, eID)
			gmsh.model.geo.addLine(cA[k] + offL, cA[k+1]+ offL, eID)
			eID += 1
	for k in range(len(cA)-1):
		print(" ADD LINE ", cA[k], cA[k] + 4, eID)
		gmsh.model.geo.addLine(cA[k], cA[k] + 4, eID)
		eID += 1

	## THE CURVES
	print(" CURVE LOOP 1 ",[1,2,3,4],  cID)  # bottom
	gmsh.model.geo.addCurveLoop([1,2,3,4],  cID)  # bottom
	print(" CURVE LOOP 2", [5,6,7,8],  cID+1)  # top
	gmsh.model.geo.addCurveLoop([5,6,7,8],  cID+1)  # top
	print(" CURVE LOOP 3", [9,-8,-12,4],  cID+2)  # left
	gmsh.model.geo.addCurveLoop([9,-8,-12,4],  cID+2)  # left
	print(" CURVE LOOP 4", [1,10,-5, -9], cID+3)  # front
	gmsh.model.geo.addCurveLoop([1,10,-5, -9], cID+3)  # front
	print(" CURVE LOOP 5", [2,11,-6, -10], cID+4)  # right
	gmsh.model.geo.addCurveLoop([2,11,-6, -10], cID+4)  # right
	print(" CURVE LOOP 6", [-3, 11,7,-12], cID+5)  # back
	gmsh.model.geo.addCurveLoop([-3, 11,7,-12], cID+5)  # back
	cID += 6

	for j in range(6):
		print(" ADD SURFACE ", [sID + j], sID + j)
		sj = gmsh.model.geo.addPlaneSurface([sID + j], sID + j)

	print(" SURFACE LOOP ", [a + sID for a in a6])
	loop = gmsh.model.geo.addSurfaceLoop([a + sID for a in a6])
	vol = gmsh.model.geo.addVolume([loop], tag=1)
	gmsh.model.geo.synchronize()

	for j in range(6):
		print(" Adding boundary ",  [sID + j], "box_"+b_names[j+1] + "_" + str(j+1))
		gmsh.model.addPhysicalGroup(2, [sID + j], name=b_names[j+1])
	sID += 6
	gmsh.model.addPhysicalGroup(3, [vol], name="volume1", tag = vol)
	# We can then generate a 2D mesh...
	gmsh.model.mesh.generate(3)
	gmsh.model.geo.synchronize()
	gmsh.write(gmsh_mesh)
	gmsh.write(gmsh_mesh.replace(".msh", ".vtk"))
	gmsh.finalize()

if create_mfem:
	
	mesh = mfem.Mesh(gmsh_mesh)
	device = mfem.Device("cpu")
	device.Print()
	dim = mesh.Dimension()
	order = 1
	fec = mfem.H1_FECollection(order, dim)
	fespace = mfem.FiniteElementSpace(mesh, fec, dim)

	# 6. Determine the list of true (i.e. conforming) essential boundary dofs.
	#    In this example, the boundary conditions are defined by marking only
	#    boundary attribute 1 from the mesh as essential and converting it to a
	#    list of true dofs.
	ess_tdof_list = mfem.intArray()
	ess_bdr = mfem.intArray([0] * (mesh.bdr_attributes.Max()))

	# impose Dirichlet BCs in the bottom box to fix objects in space
	ess_bdr[0] = 1

	fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list)

	# 7. Set up the linear form b(.) which corresponds to the right-hand side of
	#    the FEM linear system. In this case, b_i equals the boundary integral
	#    of f*phi_i where f represents a "pull down" force on the Neumann part
	#    of the boundary and phi_i are the basis functions in the finite element
	#    fespace. The force is defined by the VectorArrayCoefficient object f,
	#    which is a vector of Coefficient objects. The fact that f is non-zero
	#    on boundary attribute 2 is indicated by the use of piece-wise constants
	#    coefficient for its last component.
	f = mfem.VectorArrayCoefficient(dim)
	for i in range(dim - 1):
		f.Set(i, mfem.ConstantCoefficient(0.0))

	# this is a boundary vector. we are applying a vertical load (?)
	pull_force = mfem.Vector([0] * mesh.bdr_attributes.Max())
	pull_force[1] = -2e5

	f.Set(dim - 1, mfem.PWConstCoefficient(pull_force))
	b = mfem.LinearForm(fespace)
	b.AddBoundaryIntegrator(mfem.VectorBoundaryLFIntegrator(f))
	b.Assemble()
	# 8. Define the solution vector x as a finite element grid function
	#    corresponding to fespace. Initialize x with initial guess of zero,
	#    which satisfies the boundary conditions.
	x = mfem.GridFunction(fespace)
	x.Assign(0.0)

	# 9. Set up the material properties: top box will be rigid, bottom will be very elastic

	lamb = mfem.Vector(mesh.attributes.Max())
	mu = mfem.Vector(mesh.attributes.Max())
	# CREATE AN ELASTIC MATERIAL
	E_0 = 5e5
	nu_0 = 0.4
	lamb[0] = E_0 * nu_0 / ((1 + nu_0) * (1 - 2 * nu_0))
	mu[0] = E_0 / (2 * (1 + nu_0))
	lambda_func = mfem.PWConstCoefficient(lamb)
	mu_func = mfem.PWConstCoefficient(mu)
	a = mfem.BilinearForm(fespace)
	a.AddDomainIntegrator(mfem.ElasticityIntegrator(lambda_func, mu_func))

	# 10. Assemble the bilinear form and the corresponding linear system,
	#     applying any necessary transformations such as: eliminating boundary
	#     conditions, applying conforming constraints for non-conforming AMR,
	#     static condensation, etc.
	print('matrix...')
	a.Assemble()

	A = mfem.OperatorPtr()
	B = mfem.Vector()
	X = mfem.Vector()
	a.FormLinearSystem(ess_tdof_list, x, b, A, X, B)
	print('...done')  # Here, original version calls heigth, which is not defined in the header...!?

	# 10. Solve
	AA = mfem.OperatorHandle2SparseMatrix(A)
	M = mfem.GSSmoother(AA)
	mfem.PCG(AA, M, B, X, 1, 5000, 1e-8, 1e-8)

	# 11. Recover the solution as a finite element grid function.
	a.RecoverFEMSolution(X, b, x)

	# 13. For non-NURBS meshes, make the mesh curved based on the finite element
	#     space. This means that we define the mesh elements through a fespace
	#     based transformation of the reference element. This allows us to save
	#     the displaced mesh as a curved mesh when using high-order finite
	#     element displacement field. We assume that the initial mesh (read from
	#     the file) is not higher order curved mesh compared to the chosen FE
	#     space.
	if not mesh.NURBSext:
		mesh.SetNodalFESpace(fespace)

	# 14. Save the displaced mesh and the inverted solution (which gives the
	#     backward displacements to the original grid). This output can be
	#     viewed later using GLVis: "glvis -m mesh.mesh -g data.gf".

	nodes = mesh.GetNodes()
	nodes += x
	x.Neg()

	output = io.StringIO()
	output.precision = 8
	ref = 1
	mesh.PrintVTK(output)

	x.SaveVTK(output, 'U', ref)
	fid = open("converted_mesh.vtk", "w")
	fid.write(output.getvalue())
	fid.close()

@v-dobrev
Copy link
Member

Currently, there is no version of GridFunction::SaveVTK that works with the method Mesh::PrintVTK(std::ostream &). It is probably not too hard to add one, assuming that the GridFunction uses the same linear/quadratic order basis as the mesh.

Let us know if you want to try implementing GridFunction::SaveVTK(std::ostream &out, const std::string &field_name) that matches Mesh::PrintVTK(std::ostream &) and have any questions.

@XuliaDS
Copy link
Author

XuliaDS commented Jan 15, 2024

yep, I think it would be very useful. OK, I am currently using the pyMFEM wrap so I should implement it i C++ and then add the python call, yes ?

UPDATE. This was rather easy to write as it was already answered here: mfem/mfem#1398 (comment)

I am sorry. I didn't find that issue when I was browsing through MFEM. I rewrote the function in python:

def grid_function_save_vtk(gf: mfem.GridFunction, os: io, field_name: str):

	fes = gf.FESpace()
	mesh = fes.GetMesh()
	vec_dim = gf.VectorDim()
	assert gf.Size()/vec_dim == mesh.GetNV()
	os.write("POINT_DATA "+str(fes.GetNDofs())+"\n")
	if vec_dim == 1:
		print('Scalar field', field_name)
		os.write("SCALARS "+field_name+ " double 1\n LOOKUP_TABLE default\n")
		for i in range(fes.GetNDofs()):
			os.write(gf[i]+"\n")
	elif (vec_dim == 2 or vec_dim == 3) and mesh.SpaceDimension() > 1:
		print('Vector field', field_name)
		os.write("VECTORS "+field_name+" double\n")
		vdofs = mfem.intArray(vec_dim)
		for i in range(fes.GetNDofs()):
			vdofs.SetSize(1)
			vdofs[0] = i
			fes.DofsToVDofs(vdofs)
			os.write(str(gf[vdofs[0]]) +' '+str(gf[vdofs[1]]) +' ')
			if vec_dim == 2:
				os.write("0.0\n")
			else:
				os.write(str(gf[vdofs[2]])+"\n")
	else:
		pass

#       // other data: save the components as separate scalars
#       for (int vd = 0; vd < vec_dim; vd++)
#       {
#          out << "SCALARS " << field_name << vd << " double 1\n"
#              << "LOOKUP_TABLE default\n";
#          Array<int> vdofs(vec_dim);
#          for (int i=0; i < fes->GetNDofs(); ++i)
#          {
#             vdofs.SetSize(1);
#             vdofs[0] = i;
#             fes->DofsToVDofs(vdofs);
#             out << gf[vdofs[vd]] << '\n';
#          }
#       }
#    }

# }

and adding that to the previous script:

	output = io.StringIO()
	output.precision = 8
	ref = 1
	mesh.PrintVTK(output)
	grid_function_save_vtk(x, output, 'U')
	fid = open(os.path.join("vtkdata.vtk"), "w")
	fid.write(output.getvalue())

gives this:

image

which is exactly what I wanted :)

I will try to do the same for CELL_DATA to have it complete.

Many thanks !

Xulia

@v-dobrev
Copy link
Member

Good find @XuliaDS! We should probably add that function, SaveLinearGridFunctionVTK from mfem/mfem#1398 (comment), to the library as GridFunction::SaveVTK(std::ostream &out, const std::string &field_name). @pazner, what do you think? Also, maybe we can support quadratic H1 spaces too, since Mesh::PrintVTK supports quadratic meshes?

@XuliaDS
Copy link
Author

XuliaDS commented Jan 19, 2024

that sounds great.

FYI: one thing that I noticed is that when you visualize the mesh after running mfem, the internal patches are not shown.

Original mesh:

image

The mfem solution:

image

without full opacity:

image

In my case, what I ended up doing was then run a paraview routine that selects each cells by material. I separate those cells and end up with three volumes. Then I append them and I can see the contact regions:

image

but maybe there is better way of doing this internally already.

Cheers

Xulia

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants