vtk output for solutions

[XuliaDS]

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.

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

Thank you !

[v-dobrev]

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]

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]

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]

@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]

@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]

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]

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)

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]

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]

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]

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

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.

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]

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]

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: https://github.com/mfem/mfem/issues/1398#top

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:

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]

Good find @XuliaDS! We should probably add that function, SaveLinearGridFunctionVTK from https://github.com/mfem/mfem/issues/1398#issuecomment-610031344, 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]

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:

The mfem solution:

without full opacity:

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:

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

Cheers

Xulia