Unable to read .msh mesh file which contains only rectangular elements

mesh
22

I cannot reproduce this with replacing:

with " "

l_out = 1.0;
l_in = 0.1;
// Point of the domain
Point(1) = {0, 0, 0, l_out};
Point(2) = {10, 0, 0, l_out};
Point(3) = {10, 5-0.1, 0, l_in};
Point(4) = {10, 5+0.1, 0, l_in};
Point(5) = {10, 10, 0, l_out};
Point(6) = {0, 10, 0, l_out};
Point(7) = {0, 5+0.1, 0, l_in};
Point(8) = {0, 5-0.1, 0, l_in};
// Lines between points
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 1};
Line(9) = {8, 3};
Line(10) = {7, 4};
// Curve Loops
Curve Loop(1) = {1, 2, -9, 8}; // Line Loop
Curve Loop(2) = {4, 5, 6, 10};
Curve Loop(3) = {9, 3, -10, 7};
// Plane Surface
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
// Physical groups
// Lines
Physical Curve("L_bottom", 1) = {8, 1, 2};
Physical Curve("L_top", 2) = {4, 5, 6};
Physical Curve("L_left", 3) = {7};
Physical Curve("L_right", 4) = {3};
// Surfaces
Physical Surface("Bottom") = {1};
Physical Surface("Top") = {2};
Physical Surface("Middle") = {3};
23

Hello ! Thank you for your answer. I’m not sure I understand, do you mean that you need another file to reproduce the error?

24

I cannot reproduce the error you are showing when using the code I posted in the previous reply. The mesh is converted to xdmf without any errors

25

Unfortunately, when copying the code you posted the error remains exactly the same. Could this be due to the generation of the .msh file from gmsh? What I do is first import the .geo file into gmsh, then generate the 2d mesh (mesh->define->2d), then smooth it (mesh->define->smooth 2d). Finally, I export the file as an .msh file. For the export, I use the Version 4 ASCII format and save all elements. Am I making a mistake somewhere?

26

You should not use SaveAll.
Note that you can add smoothing options and SaveAll 0 to your geo file:

Mesh.Smoothing = 100;
Mesh.SaveAll = 0;

You can set the smoothing parameter to any integer.

With these two additional lines in your geo file, you can simply create the mesh by calling gmsh -2 mesh.geo

1 Like
27

Thank you very much for your reply. It worked fine!

I have a second question regarding the application of boundary conditions. Once I get the cell_tags and facet_tags with the following code (see below), how can I apply boundary conditions?

# Read the mesh
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    cell_tags = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

with XDMFFile(MPI.COMM_WORLD, "facet_mesh.xdmf", "r") as xdmf:
    facet_tags = xdmf.read_meshtags(mesh, name="Grid")

To take a concrete example, in my case, I have a square mesh in which each side is a physical curve: L_bottom, L_top, L_left and L_right. There is also a physical surface which represents the whole domain (its surface). If for example I want to apply a Dirichlet boundary condition on L_left and a different one on all other sides, how can I get these physical curves to apply these conditions ? More generally, I’m having a bit of trouble understanding what these MeshTags_int32 objects are.
Thanks again for your help!

28

See for instance: Combining Dirichlet and Neumann conditions — FEniCSx tutorial or any of the other tutorials regarding boundary conditions in my tutorial

29

Thank you! So if I understood correctly, it is possible to get the numbers of the nodes that are on the border by accessing the facet_tags.indices attribute, and it is possible to see to which physical group it belongs by accessing the facet_tags.values attribute. However, the facet_tags.values attribute returns integers (for example between 1 and 4 if I have defined 4 physical curves for each of the edges of the square). Is there a way to select the nodes of a certain group directly by the name that was given to the group in the .geo file?

For the moment I do facet_tags.indices[facet_tags.values == 1] to get the nodes of group 1 but I don’t know if this group 1 corresponds to the physical curve L_top, L_bottom, L_left or L_right.

30

You should define the surfaces in gmsh with an integer tag as well as with a string. We do not support reading in tags defined through strings, as strings can be of arbitrary size, and comparing them is a more costly operation than comparing integers. Note that you can open the mesh (.msh) in GMSH to see what values you have tagged each surface with. You can also do this by opening the .xdmf file in for instance Paraview.

1 Like