Getting salome mesh to FEniCSx python

mesh
1

Hello,
first of all, my excuses if i am not clear enought please bear with me as It is my first post.

I want to use FEniCSx to do some FEM simulations, I do not use Gmsh but Salome for my geometry and mesh generation.
I installed FEniCSx using:

sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt update
sudo apt install fenicsx

and now in python i can import dolfinx

The first thing i would like to do is be able to import my mesh generated in salome:

In salome, similar to gmsh, i can create dedicated groups for the boundary conditions and other, 1D, 2D and 3D groups (where they will contain, edges, faces, and volumes/cells of the mesh).
Nativally I can export the mesh in the following formats: GMF Ascii, MED 4.2, IDEAS (unv),cgns.
I could also use the integrated meshio to export to .msh (nevertheless I am less keen to use the integrated meshio, as it has proved in the past to be buggy, such as it has the exporting on format xdmf which from what I found is the ideal format for fenicsx does not work.).

Reason for this current post,
I am looking to bring the mesh created (with the different groups) from salome to a fenicsx simulation. I tried using (in a python file test.py outside of salome):

import meshio
mesh = meshio.read(
    '/home/franco/Downloads/dummyFiles/test/hexa.med'
)
mesh.write(
    '/home/franco/Downloads/dummyFiles/test/hexa.xdmf' )

Then run it by using python3 test.py, this proven to ‘write’ the mesh in the desired xdmf format, nevertheless, i am having some trouble still:
The groups I have created in salome disappeared, there is only a cell_tags i can see in paraview, that are negative from -12 to 0, where each negative int value looks like it corresponds to one of the groups that i had in salome, from the -12 until the -6 but the -5 to -1 are empty and furthermore the 0 I dont finish of recognizing what corresponds to as it does not completely correspond to anything I have a tag on salome nor the submeshes.

my questions and what i would require help/recommendations would be:

  1. with the formats i can export to, what would be a good workflow where I could keep the boundary names? to use in the simulation after? for example, I can see that the UNV format keeps the names of the groups stored in the file, I can not say it for the MED nor cgns (as they are write in binary), the GMF does not save anything related to the groups,
  2. once I get the mesh in the correct format with all the information (point 1.) how can I import it to fenicxs to run a simple simulation over it? I did not find any tutorial with mesh importation + definition of a simple boundary conditions.

I am sorry if it is not clear enough, I do not have experience with fenicxs (nor gmsh), but have some with OpenFOAM and a lot with salome. what I would like to succed right now is to do the complete pre processing with salome for fenicxs as I did for OpenFOAM, that once the workflow was developed i could simply ‘plug’ it to it and get my simulation to be able to run (obviously after one needs to define the boundary conditions and physics for the solver).

thanks in advance,

2

I usually do more than what you have done here when working with meshio.
An example of this (using gmsh) can be found in many posts in the forum, such as):

which points to my tutorial and is similar to:

Note that I’ve been rewriting adios4dolfinx to support all kinds of IO backends, including:

  • pyvista (meshio)
  • VTKHDF
  • XDMF

in: General backend support by jorgensd · Pull Request #193 · jorgensd/adios4dolfinx · GitHub

I could help you read the .med file if you provide an example.
For instance I made a bespoke writer for the exodus file format at: GitHub - jorgensd/mesh_converter: A mesh converter from EXODUS 2 to XDMF

There are many demos that does this, for instance at:

3

god, awsome answer! thanks a lot!

of course, for the .med file here you have two examples, one 2D and one 3D.
here you have a mesh in 3D (with the boundary groups created FileSender)
and here a 2D version (FileSender)
keep in mind that the med file is the default output of salome so it could be interesting for fenicxs to have it, as salome inside of it has gmsh mesher and other meshers.
I will check all the links you gave me, thanks again!

4

Hello again,
I have some questions in regards to the first link you send me:
in the meshio.write() section it looks like one needs to:

  1. write the mesh meshio.write("mesh.xdmf", triangle_mesh, compression=None)
  2. write a separate file for the groups? (instead of being inside of the mesh file?) meshio.write("mt.xdmf", line_mesh, compression=None)
    so if i have a mesh of a cube (therefore 3D) and i have 6 different boundaries conditions (for each face of the cube). I should write the mesh of the cube, tetra/hexa or whatever, and also write 6 different files for each boundary? so in total I would need to create 7 files, and the ones of the boundaries would have repeated information? as the cube would have this faces inside?
    i am understanding it wrongly?
    2.b. if we have different groups in 3D, eg. we have a mesh with half of it for one domain and the other half for the other domain (eg., air and solid) I need to: 1. write a volumetric mesh with all the elements, 2. write another file ONLY containing the elements of the first domain and anotherone for the other domain? so basically doubling the mesh (in files I mean)?
  3. also, the function I am seeing is for a monotype mesh (ie., only one kind of element), but how can we do for mix element meshes? eg. tetrahedral+pyramids+hexahedral meshes?
5

I started on a MED-reader at:

it can currently read in mixed-cell-type meshes to DOLFINx (I tested on your two med files).
Note that it cannot read in the cell or facet markers (as I would have to change some internal logic in DOLFINx to make this work for mixed-cell-type grids.

If you have a single cell type grid (pure quad, pure triangle, pure tet or pure hex), I can extend the reader to read in the different physical sub-groups.

No, for meshio’s/DOLFINx XDMF-format to work, we have to write each cell type to a separate file. I.e. write your cells to one file (tetra or hex), then write all facet markers to another file (facet makers), (one for edges and one for ridges and one for peaks, all these are optional).

No, one would write all cell markers within the same file, it is only the cell-type that matters.

Currently there is limited supported for mixed element meshes.
There is a demo at:
Poisson equation — DOLFINx 0.10.0.post5 documentation that uses hexahedral and prism elements (but the interface isn’t super nice).

Currently, one can read in such meshes from file with dolfinx.io.vtkhdf.read_mesh (dolfinx.io.vtkhdf — DOLFINx 0.10.0.post5 documentation) which supports the VTKHDF format (VTKHDF File Format - VTK documentation)

2 Likes
6

thanks for your answer and help, i was just about to beging in writing my own mesh parser inside salome, by any chance could i ask for a small example of a mesh in the ‘perfect’ format for DOLFINx? 2D and 3D (with all the files i mean facet and others). if it is not much to bother.

it can currently read in mixed-cell-type meshes to DOLFINx (I tested on your two med files).
Note that it cannot read in the cell or facet markers (as I would have to change some internal logic in DOLFINx to make this work for mixed-cell-type grids.

the the cell or facet markers would be to ‘indicate’ the physical groups, no?

If you have a single cell type grid (pure quad, pure triangle, pure tet or pure hex), I can extend the reader to read in the different physical sub-groups.

here you have 2 meshes, one quad with the 1D groups for the boundary conditions, and the tetra mesh with 2D groups for the boundary conditions and 3D groups for the two regions:

No, for meshio’s/DOLFINx XDMF-format to work, we have to write each cell type to a separate file. I.e. write your cells to one file (tetra or hex), then write all facet markers to another file (facet makers), (one for edges and one for ridges and one for peaks, all these are optional).

again, the markers are the indications for the boundary condition and regions, no?

Currently there is limited supported for mixed element meshes.

yes i discovered this yesterday, i thought it was not the case, reason for the meshes i shared originally, i will stay in pure meshes for the moment (as i dont require mixed meshes for the moment)

7

Yes. In general, DOLFINx supports storing all of these in a single XDMF file, but this is not the case when using meshio, as it doesn’t support that specific XDMF specification.

I’ve made a reader that first reads the mesh data directly into DOLFINx (including cells and facet markers) and stores it in the XDMFFile format, that is currently the preferred reader of meshes and tags (will hopefully be replaced by the VTKHDF format soon).
The reader is just a few more lines to the code I initially shared. The new code is available at the same link:

The code also shows how to read in the produced .xdmf file in DOLFINx (which works in parallel).

I will probably include this reader inside my tool adios4dolfinx (soon io4dolfnx) as part of: General backend support by jorgensd · Pull Request #193 · jorgensd/adios4dolfinx · GitHub

In general DOLFINx supports any format, as it can take in raw data (as shown in the aforementioned link). The two default formats we support is .msh and .xdmf, which in general can store all the mesh data (including domain markers) within a single file.

1 Like
8

Note that the latest push to the gist makes it work with arbitrarily many MPI processors. The cells and facet (markers) are read in as distributed datasets, while the nodes are read in on the first process and then distributed with DOLFINx (as they need a re-numbering according to the .med “NUM” key).

1 Like
9

Again, thanks for your help,
might be important to point, not sure if it is in the MED format also, but salome uses a different numbering for the nodes of cells in comparation to vtk (which i assume from chatgpt) that is the ordering of xdmf format.

I was writing my own mesh writer to xdmf file directly, but I am having some troubles, would you mind give me a feed back?


def exportXdmf(Mesh,path,groupsSurface=[]):
    defaultPrecisionInt=4
    defaultPrecisionFloat=8
    nodesIDs = Mesh.GetNodesId()
    facesIDs = Mesh.GetElementsByType(SMESH.FACE)
    # cellsIDs = Mesh.GetElementsByType(SMESH.VOLUME)
    offsetID=1
    #groupsSurface
    groupsIDsSurf=[g.GetIDs() for g in groupsSurface]
    groupsNamesSurf=[g.GetName() for g in groupsSurface]
    groupsIDsSurfExp=set(sum(groupsIDsSurf,[]))
    #2D section
    typesOfFaces=[SMESH.Geom_TRIANGLE, SMESH.Geom_QUADRANGLE]
    namesFaces=["Triangle","Quadrilateral"]
    orderTri=[0,1,2] #correct order
    orderQuad=[0,1,2,3] #correct order
    ordersFaces=[orderTri,orderQuad]
    numberOfElemsFaces=[Mesh.NbTriangles(), Mesh.NbQuadrangles()]
    #3D section
    typesOfCells=[SMESH.Geom_TETRA, SMESH.Geom_PYRAMID, SMESH.Geom_PENTA, SMESH.Geom_HEXA]
    namesCells=["Tetrahedron","Pyramid","Wedge","Hexahedron"]
    orderTetra=[0, 2, 1, 3] #[0,1,2,3]
    orderPyr=[0, 3, 2, 1, 4] #[0,1,2,3,4]
    orderWed=[0,1,2,3,4,5] #correct order
    orderHexas=[0, 3, 2, 1, 4, 7, 6, 5] #[1,5,4,0,2,6,7,3]
    ordersCells=[orderTetra,orderPyr,orderWed,orderHexas]
    numberOfElemsCells=[Mesh.NbTetras(), Mesh.NbPyramids(), Mesh.NbPrisms(), Mesh.NbHexas()]
    lines=[f"<?xml version=\"1.0\" ?>"]
    lines+=[f"<Xdmf Version=\"3.0\">"]
    lines+=[f"  <Domain>"]
    lines+=[f""]
    lines+=[f"    <!-- ===================================================== -->"]
    lines+=[f"    <!--                       Main Mesh                       -->"]
    lines+=[f"    <!-- ===================================================== -->"]
    lines+=[f""]
    if sum(numberOfElemsCells)!=0:
        lines+=[f"    <Grid Name=\"volume_mesh\" GridType=\"Uniform\">"]
        lines+=[f""]
        for nType, name in enumerate(namesCells):
            if numberOfElemsCells[nType]:
                cellsIDs=Mesh.GetIdsFromFilter(smesh.GetFilter(SMESH.VOLUME, SMESH.FT_ElemGeomType, typesOfCells[nType]))
                lines+=[f"      <Topology TopologyType=\"{name}\" NumberOfElements=\"{numberOfElemsCells[nType]}\">"]
                lines+=[f"        <DataItem Format=\"XML\""]
                lines+=[f"                  Dimensions=\"{numberOfElemsCells[nType]} {len(ordersCells[nType])}\""]
                lines+=[f"                  NumberType=\"Int\""]
                lines+=[f"                  Precision=\"{defaultPrecisionInt}\">"]
                for cellID in cellsIDs:
                    nodesIDsOnCell=[ID-offsetID for ID in Mesh.GetElemNodes(cellID)]
                    nodesIDsOnCell=[str(nodesIDsOnCell[n]) for n in ordersCells[nType]]
                    lines+=[f"          {' '.join(nodesIDsOnCell)}"]
                lines+=[f"        </DataItem>"]
                lines+=[f"      </Topology>"]
    lines+=[f""]
    lines+=[f"      <Geometry GeometryType=\"XYZ\">"]
    lines+=[f"        <DataItem Format=\"XML\""]
    lines+=[f"                  Dimensions=\"{len(nodesIDs)} 3\""]
    lines+=[f"                  NumberType=\"Float\""]
    lines+=[f"                  Precision=\"{defaultPrecisionFloat}\">"]
    for nodeID in nodesIDs:
        xyz=[round(i,8) for i in Mesh.GetNodeXYZ(nodeID)]
        lines+=[f"          {xyz[0]} {xyz[1]} {xyz[2]}"]
    lines+=[f"        </DataItem>"]
    lines+=[f"      </Geometry>"]
    if sum(numberOfElemsFaces)!=0:
        lines+=[f""]
        lines+=[f"    </Grid>"]
        lines+=[f""]
        lines+=[f"    <!-- ===================================================== -->"]
        lines+=[f"    <!--                Boundary Surface Mesh                  -->"]
        lines+=[f"    <!-- ===================================================== -->"]
        lines+=[f""]
        lines+=[f"    <Grid Name=\"boundary_faces\" GridType=\"Uniform\">"]
        lines+=[f""]
        for nType, name in enumerate(namesFaces):
            if numberOfElemsFaces[nType]:
                cellsIDs=set(Mesh.GetIdsFromFilter(smesh.GetFilter(SMESH.FACE, SMESH.FT_ElemGeomType, typesOfFaces[nType])))
                cellsIDs=cellsIDs & groupsIDsSurfExp
                lines+=[f"      <Topology TopologyType=\"{name}\" NumberOfElements=\"{len(cellsIDs)}\">"]
                lines+=[f"        <DataItem Format=\"XML\""]
                lines+=[f"                  Dimensions=\"{len(cellsIDs)} {len(ordersFaces[nType])}\""]
                lines+=[f"                  NumberType=\"Int\""]
                lines+=[f"                  Precision=\"{defaultPrecisionInt}\">"]
                for groupIDs in groupsIDsSurf:
                    for cellID in groupIDs:
                        nodesIDsOnCell=[ID-offsetID for ID in Mesh.GetElemNodes(cellID)]
                        nodesIDsOnCell=[str(nodesIDsOnCell[n]) for n in ordersFaces[nType]]
                        lines+=[f"          {' '.join(nodesIDsOnCell)}"]
                lines+=[f"        </DataItem>"]
                lines+=[f"      </Topology>"]
    lines+=[f""]
    lines+=[f"      <!-- Reuse same geometry -->"]
    lines+=[f"      <Geometry Reference=\"/Xdmf/Domain/Grid[@Name='volume_mesh']/Geometry\"/>"]
    lines+=[f""]
    lines+=[f""]
    lines+=[f""]
    lines+=[f"      <!-- Face Tags -->"]
    lines+=[f"      <Attribute Name=\"face_tags\""]
    lines+=[f"                 AttributeType=\"Scalar\""]
    lines+=[f"                 Center=\"Cell\">"]
    lines+=[f"        <DataItem Format=\"XML\""]
    lines+=[f"                  Dimensions=\"{len(groupsIDsSurfExp)}\""]
    lines+=[f"                  NumberType=\"Int\""]
    lines+=[f"                  Precision=\"{defaultPrecisionInt}\">"]
    lines+=[f""]
    for nG,groupIDs in enumerate(groupsIDsSurf):
        lines+=['          '+" ".join([str(nG+1) for _ in groupIDs])]
    lines+=[f""]
    lines+=[f"        </DataItem>"]
    lines+=[f"      </Attribute>"]
    lines+=[f""]
    lines+=[f"      <!-- Tag meaning metadata -->"]
    lines+=[f"      <Information Name=\"tagMeaning\">"]
    for nG,name in enumerate(groupsNamesSurf):
        lines+=[f"        {nG+1} = {name}"]
    lines+=[f"      </Information>"]
    lines+=[f"    </Grid>"]
    lines+=[f"  </Domain>"]
    lines+=[f"</Xdmf>"]
    with open(path, 'w', encoding='utf-8') as f:
        f.writelines(line + '\n' for line in lines)
    return

it is currently creating the mesh, but i see in paraview some issues that i think it is due to the ordering of the nodes of the elements maybe? i see overlapping surfaces like this:

image
image532×446 6.23 KB

salome ordering is different to the one of vtk (that is the way xdmf stores it no?)
https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
Nodal connectivity of elements — Mesh 9.15.0 documentation

here is an example of the output mesh i am getting FileSender
in paraview if i use cell size to check the volume it looks like it is positive for the cells as expected. (notice that if i use the salome default ordering for the tetra, hexa and pyr which i modify in the function, i would get negative volumes)

10

For what cell types? All the meshes I tested that you have given me works nicely with DOLFINx.

I don’t really have the bandwidth to make comments on alternative solutions to what I’ve already shown you.
Especially as I have no clue what class the Mesh in:

is.

A general note is that one should rather use ElementTree and h5py to write XDMFFile’s as for instance done in: scifem/src/scifem/xdmf.py at b06118b09df7923fa72190795829b68cf5bc63a2 · scientificcomputing/scifem · GitHub

In general, I never trust Paraview to compute cell volumes, as their specific ordering to compute volumes relies on the possibility of globally re-orienting every cell.

11

the hexahedrons/pyramids/tetrahedrals needs a re ordering if one wants to follow the vtk formating at least.

I don’t really have the bandwidth to make comments on alternative solutions to what I’ve already shown you.

thanks in any case, I just wanted to not go thought an intermediate format, and write it directly in the correct one. Mesh class is from salome itself, I have enought knowledge on it to be able to get the information i require but my issue is with the xdmf file format that i could not find a source (like the one i posted for vtk in my previous comment).

A general note is that one should rather use ElementTree and h5py to write XDMFFile’s as for instance done in: scifem/src/scifem/xdmf.py at b06118b09df7923fa72190795829b68cf5bc63a2 · scientificcomputing/scifem · GitHub

I will check this out, thanks.

12

You are not going through an intermediate format with what I am proposing in my Github Gist as it reads the data directly into DOLFINx.

The hexahedron seems to be ok to me in your mesh:

Have you tried using the code i propose for the meshes you sent and looked at the results?

1 Like
13

yes of course! but i began implementing my code before you posted yours, and I am not going to be the one using it, i am in charge of the meshing part (at least for the moment) so i need to develope this pipeline, and the thing is when you posted yours, i went though it and i did not understood it as you are using a lot of python tools that I do not know, and i am afraid of plugging this in the pipeline and then people comming and saying ‘it does not work i dont know why’ or i can not use it as intended and then not having an answer for them. as when you posted i was almost at the end of finishing my script i continue yesterday to finish it.

14

Have you tried using the code i propose for the meshes you sent and looked at the results?

also, I just exported the salome Mesh object into med, and used your script from here and when i try to open the test_med.xdmf file with paraview, it directly crash paraview.
here you have the med file i tested with. FileSender

and if i may, what is the command you used to read the xdmf file I posted?
when i do (following the links you post it in the begining)

with XDMFFile(MPI.COMM_WORLD, "/home/franco/Downloads/dummyFiles/test/tq.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")

i am getting:

Traceback (most recent call last):
  File "/home/franco/Downloads/dummyFiles/test/Mesh-1.py", line 25, in <module>
    with XDMFFile(MPI.COMM_WORLD, "/home/franco/Downloads/dummyFiles/test/tq.xdmf", "r") as xdmf:
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Unable to open HDF5 file. File /home/franco/Downloads/dummyFiles/test/tq.h5 does not exist.

I mean the mesh i post it it does not have the cell tags in another file (tq.h5) as it is classically done, reason why fenicxs is not linking it, how did you conturnate this?
thanks in advance

15

You should use the extract_blocks filter in paraview to look at each individual component. Here is how the tq.xdmf file you sent looks like:

image
image1920×1052 480 KB

You have written a plain-text XDMF file, rather than a XDMF file with pointers to a HDF5 encoded file, and should therefore specify ASCII encoding in the initiation of XDMFFile, see: dolfinx.io — DOLFINx 0.10.0.post5 documentation
The default encoding is HDF5, as it is a scalable format that can easily be read in parallel.
This is similar to what the .med format does, as .med is just a specification of a layout within the HDF5-encoded files.

What XDMF3 Reader did you choose when opening the file in Paraview?
One should select the XDMF3ReaderT, as shown in the screenshot below.

image
image508×530 25.7 KB

Converting your tq_.med file gives the result below (not that I clipped the facet tags to show the interior tag.

image
image1920×1052 269 KB

1 Like
16

You should use the extract_blocks filter in paraview to look at each individual component. Here is how the tq.xdmf file you sent looks like:

yes i know the extract block filter, but in comparation to the xdmf files i have open as examples (even yours) do not look like there is an overlapping surfaces, it looks like it is an issue from normals but i could not find how to solve it for the moment.

You have written a plain-text XDMF file, rather than a XDMF file with pointers to a HDF5 encoded file, and should therefore specify ASCII encoding in the initiation of XDMFFile , see: dolfinx.io — DOLFINx 0.10.0.post5 documentation

great! thanks i will check it!

What XDMF3 Reader did you choose when opening the file in Paraview?
One should select the XDMF3ReaderT , as shown in the screenshot below.

i was using the xdmf reader directly not the XDMF3ReaderT, strange as the previous meshes where not crashing.

thanks a lot for your help,
regards,

17

Where above do you se overlap? I don’t see any obvious overlap there.

18

Where above do you se overlap? I don’t see any obvious overlap there.

you misunderstood me, i meant with the mesh I posted myself, tq.xdmf not the one produced by the med reader script.

image
image669×483 16.3 KB

19

But this is before you extract blocks? If I use the normal XMDFReader (not the T-reader) on your tq.xdmf I can extract each block to see that there is no overlap.

20

No i did not, but in other files (that i did not make myself) i got as example without extracting the blocks it did not showed. and as i have writen in the past some mesh parsers i have seen this kind of ‘visual effect’ when the normals were wrong (ie ordering of the nodes on faces or cells)