Using exported XDMF mesh file from FreeCAD on FEniCS

[tschnk]

Hi all,

I am running the current (stable) Version of FreeCAD 0.18 (Rev. 16146 Git) as an Appimage on Linux Mint 20 which works great so far and created a simple box as geometry and meshed it with the included gmsh mesher. After that I’ve created mesh groups to mark one face of the box and another mesh group to mark the face on the opposite. I’ve also marked the whole solid and so I got three groups, related to the whole mesh. They all got a own defined label and I’ve set “Use Label” to “true” on each group.

Now my need is to find any workflow so I could use them on a FEniCS Python script for integrating over subdomains, defined in FreeCAD, to implement subdomain data like:

dV = Measure(‘dx’, domain=mesh, subdomain_data=volume subdomains, here just “solid”, metadata={‘quadrature_degree’: 2})
dA = Measure(‘ds’, domain=mesh, subdomain_data=facet subdomains, here “rear” and “front”, metadata={‘quadrature_degree’: 2})

Finally a test would be to integrate over any of these subdomains with the “assemble” function of FEniCS. So that:

assemble(1*dV(solid_sub_domain_from_freecad))

will output the volume of the box (10 mm x 10 mm x 10 mm = 1000 mm^3).

(Previously I used to import a .msh to xml converted mesh, added it with "mesh = Mesh(‘mesh.xml’) and got the Mesh Functions by using MeshFunction(‘size_t’, mesh, ‘mesh_physical_region.xml’) for the volumes and MeshFunction(‘size_t’, mesh, ‘mesh_facet_region.xml’) for the facet regions respectively.)

Is there any user who have experience in that workflow to use that XDMF mesh export file(s) efficiently on FEniCS, including subdomain data defined in FreeCAD and wants to share a minimal script to do that?

Best regards,
Tobias

[dokken]

There are many post on this throughout the forum (simply search for Measure or XDMFFile) and you Get results such as:

Gmsh mesh volume (and surface) is worth 0

I am facing another mesh problem, which I used to solve by using an old mesh format with gmsh (msh2), but this time this won’t fix the problem. I have the following mesh.geo: SetFactory("OpenCASCADE"); Box(1) = {0, 0, 0, 1e-2, 5e-3, 1e-4}; //+ Box(2) = {2.5e-3, 2.5e-3, 1e-4, 5e-3, 1e-4, 1e-4}; //+ BooleanFragments{ Volume{1}; Delete; }{ Volume{2}; Delete; } //+ Physical Volume("the_box") = {1}; //+ Physical Volume("thin box") = {2}; //+ Physical Surface("end_left") = {7}; //+ Physical Surface(…

and

Different results for the same mesh created differently (Hyperelasticity)

I claim that the meshes are not the same, see attached picture: [cube_comparison] } I would also recommend you not to use xml to load meshes, as it has been deprecated for a while. The suggested way of converting a mesh from gmsh is by using meshio. Also, to make sure that your visualization is P2 data is saved, use xdmf.write_checkpoint. Here is a script doing all of the suggested actions: from dolfin import * import meshio import numpy as np class Left(SubDomain): def inside(self, x,…

[tschnk]

First I want to thank you for your quick reply!

The main question is if it is possible to only use the .xdmf meshes generated/exported by FreeCAD without the use of converting from msh to xdmf by meshio.

[dokken]

What happens if you try to import it (using. XDMFFile.read as shown in the links).

The easiest would be to test the approaches suggested, and then report a minimal working example with an ascii based xdmf file that you can share with you error message.

There are different ways of writing xdmf, and not all of them are suitable for parallel computing.

[tschnk]

So I’ve exported the mesh of a 10 mm x 10 mm x 10 mm box (meshed with gmsh within FreeCAD). As FreeCAD asked me to give self defined IDs, I have chosen individual IDs (2,4,8) for “Face”, “Solid” and “Top” respectively:

(FreeCAD users must know to use first order meshes, see Using XDMF Mesh Export File(s) for FEniCS - FreeCAD Forum)

Because of this post:
https://fenicsproject.discourse.group/t/gmsh-mesh-volume-and-surface-is-worth-0/4944

I thought that three meshes would be neccessary to have a successful workaround to handle subdomain data. Previously I used the “main mesh” xml file, a seperate file for “physical_regions” and another seperate file for “facet_regions”.

After exporting in FreeCAD I got one XDMF Mesh file I will show in this post below. A test by using the “assemble” function gives accurate values:

from dolfin import *

mesh = Mesh()
with XDMFFile("box_test.xdmf") as infile:
    infile.read(mesh)

physical_reg = MeshValueCollection('size_t', mesh, mesh.topology().dim())
with XDMFFile("box_test.xdmf") as infile:
    infile.read(physical_reg, "Solid")
cf = cpp.mesh.MeshFunctionSizet(mesh, physical_reg)

facet_reg = MeshValueCollection('size_t', mesh, mesh.topology().dim()-1)
with XDMFFile("box_test.xdmf") as infile:
    infile.read(facet_reg, "Face")
ff = cpp.mesh.MeshFunctionSizet(mesh, facet_reg)


facet_reg2 = MeshValueCollection('size_t', mesh, mesh.topology().dim()-1)
with XDMFFile("box_test.xdmf") as infile:
    infile.read(facet_reg2, "Top")
ff2 = cpp.mesh.MeshFunctionSizet(mesh, facet_reg2)


dx = Measure("dx", domain=mesh, subdomain_data=cf, metadata={'quadrature_degree': 2})
ds_front = Measure("ds", domain=mesh, subdomain_data=ff, metadata={'quadrature_degree': 2})
ds_top = Measure("ds", domain=mesh, subdomain_data=ff2, metadata={'quadrature_degree': 2})

vol = assemble(1*dx(4))
surf = assemble(1*ds_front(2))
surf_top = assemble(1*ds_top(subdomain_id=8))

print("total volume: " + str(vol))
print("front surface: " + str(surf))
print("top surface: " + str(surf_top))

which seems to look good so far. The output gives:

total volume: 999.9999999999993
front surface: 100.00000000000003
top surface: 100.00000000000003

In the next few days I may try to extend the code to solve a simple example to solve Poisson equation to give a workaround for interested who want to do the same thing someday.

So if you don’t see anything “strange” in the XDMF file, everything should be fine so far!
I will be in touch if there are any problems after all.

box_test.xdmf:

<?xml version="1.0"?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf version="3.0"><Domain><Grid Name="base_mesh" GridType="Uniform"><Topology TopologyType="tetrahedron" NumberOfElements="24" NodesPerElement="4"><DataItem NumberType="UInt" Dimensions="24 4" Format="XML">9 11 10 12
8 13 11 10
11 10 13 9
8 10 11 12
1 0 8 10
0 2 8 13
10 0 13 4
3 11 8 2
1 8 3 12
11 13 2 6
4 13 9 6
6 11 9 7
10 4 9 5
11 7 3 12
12 9 7 5
12 1 10 5
13 8 0 10
8 2 11 13
10 8 1 12
13 9 10 4
11 9 13 6
11 3 8 12
7 11 9 12
12 10 9 5</DataItem></Topology><Geometry GeometryType="XYZ"><DataItem Dimensions="14 3" Format="XML">0.000000 0.000000 10.000000
0.000000 0.000000 0.000000
0.000000 10.000000 10.000000
0.000000 10.000000 0.000000
10.000000 0.000000 10.000000
10.000000 0.000000 0.000000
10.000000 10.000000 10.000000
10.000000 10.000000 0.000000
0.000000 5.000000 5.000000
10.000000 5.000000 5.000000
5.000000 0.000000 5.000000
5.000000 10.000000 5.000000
5.000000 5.000000 0.000000
5.000000 5.000000 10.000000</DataItem></Geometry></Grid><Grid Name="Face_mesh" GridType="Uniform"><Topology TopologyType="triangle" NumberOfElements="24" NodesPerElement="3"><DataItem NumberType="UInt" Dimensions="24 3" Format="XML">0 1 10
4 0 10
1 5 10
5 4 10
0 13 2
4 13 0
2 13 6
6 13 4
1 0 8
0 2 8
3 1 8
2 3 8
5 9 4
4 9 6
7 9 5
6 9 7
2 11 3
6 11 2
3 11 7
7 11 6
1 3 12
5 1 12
3 7 12
7 5 12</DataItem></Topology><Geometry Reference="XML">/Xdmf/Domain/Grid/Geometry</Geometry><Attribute AttributeType="Scalar" Center="Cell" Name="Face"><DataItem Dimensions="24 1" Format="XML">2
2
2
2
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0</DataItem></Attribute></Grid><Grid Name="Solid_mesh" GridType="Uniform"><Topology TopologyType="tetrahedron" NumberOfElements="24" NodesPerElement="4"><DataItem NumberType="UInt" Dimensions="24 4" Format="XML">9 11 10 12
8 13 11 10
11 10 13 9
8 10 11 12
1 0 8 10
0 2 8 13
10 0 13 4
3 11 8 2
1 8 3 12
11 13 2 6
4 13 9 6
6 11 9 7
10 4 9 5
11 7 3 12
12 9 7 5
12 1 10 5
13 8 0 10
8 2 11 13
10 8 1 12
13 9 10 4
11 9 13 6
11 3 8 12
7 11 9 12
12 10 9 5</DataItem></Topology><Geometry Reference="XML">/Xdmf/Domain/Grid/Geometry</Geometry><Attribute AttributeType="Scalar" Center="Cell" Name="Solid"><DataItem Dimensions="24 1" Format="XML">4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4</DataItem></Attribute></Grid><Grid Name="Top_mesh" GridType="Uniform"><Topology TopologyType="triangle" NumberOfElements="24" NodesPerElement="3"><DataItem NumberType="UInt" Dimensions="24 3" Format="XML">0 1 10
4 0 10
1 5 10
5 4 10
0 13 2
4 13 0
2 13 6
6 13 4
1 0 8
0 2 8
3 1 8
2 3 8
5 9 4
4 9 6
7 9 5
6 9 7
2 11 3
6 11 2
3 11 7
7 11 6
1 3 12
5 1 12
3 7 12
7 5 12</DataItem></Topology><Geometry Reference="XML">/Xdmf/Domain/Grid/Geometry</Geometry><Attribute AttributeType="Scalar" Center="Cell" Name="Top"><DataItem Dimensions="24 1" Format="XML">0
0
0
0
8
8
8
8
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0</DataItem></Attribute></Grid></Domain></Xdmf>