Extracting mesh regions for gmsh mesh

mesh
1

Hi everyone,

I’m currently trying to read a gmsh-generated mesh for a cuboid channel with 3 different types of boundaries.

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
3
2 1 "ExcitedBoundary"
2 2 "SoundProofBoundary"
2 3 "AdmittanceBoundary"
$EndPhysicalNames
$Nodes
...
$EndNodes
$Elements
1052
1 2 2 1 1 2 13 3
2 2 2 1 1 3 13 12
3 2 2 1 1 9 6 7
...
11 2 2 2 2 8 5 49
12 2 2 2 2 49 5 58
13 2 2 2 2 17 6 7
14 2 2 2 2 17 90 6
...

I can convert it using

pathname = "./DuctMesh/"
msh = meshio.read(pathname + "DuctMesh.msh")

meshio.write(pathname + "DuctMesh_xdmf.xdmf", meshio.Mesh(points=msh.points,
                                                          cells={"triangle": msh.cells["triangle"]}))

meshio.write(pathname + "DuctMesh_xdmf_Boundary.xdmf",
             meshio.Mesh(points=msh.points,
                         cells={"triangle": msh.cells["triangle"]},
                         cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))

and import it to Fenics using

mesh = df.Mesh()

with df.XDMFFile(pathname + "DuctMesh_xdmf.xdmf") as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection("size_t", mesh, 2)
print("The MeshValueCollection: ", mvc)
with df.XDMFFile(pathname + "DuctMesh_xdmf_Boundary.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)

But i don’t know how to mark my boundaries properly or how to acces them from the previous code.
When I’m using

boundary_markers = df.MeshFunction("size_t", mesh(), mesh.topology().dim() - 1)

and

class OnBoundary(df.SubDomain):  # determine whether coordinate is on boundary

    def inside(self, x, on_boundary):
        return on_boundary

i am not able to mark my boundries in the right way and accesing them properly.

neumann_excited = Bm.OnBoundary()  # define Neumann boundary
neumann_excited.mark(boundary_markers, 1)  # mark excited neumann boundary as 1
neumann_zero = Bm.OnBoundary() 
neumann_zero.mark(boundary_markers, 2)  # mark zero neumann boundary as 2
robin = OnBoundary()  # define Robin boundary
robin.mark(boundary_markers, 3)  # mark robin boundary as 3

When i’m debuggin my boundaries, each point is assigned to each of my boundaries, but i want boundary one to be neumann_excited, etc. …

I am not able to find anythin on how to solve this problem, so i’m asking here and hoping i can get an answer.

Any help is highly appreciated!

Thanks in advance!

1 Like
2

Hi,
Your boundary markers are already loaded into the facet function mf.
The line mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc) create the meshfunction that can be used in for instance a measure or Dirichlet-condition:
ds_admittance = df.Measure(subdomain=mesh, subdomain_data=mf, subdomain_id=3), or
bc = DirichletBC(FunctionSpace(mesh, "CG", 1), Constant(1), mf, 3)

1 Like
3

Wow thank your very much, that was very useful.
Do i even have to use this class definition

class OnBoundary(df.SubDomain):  # determine whether coordinate is on boundary

    def inside(self, x, on_boundary):
        return on_boundary

if i’m defining my integrals to be

ds_excited = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=1)
ds_sound_proof = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_admittance = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)

And just for understanding what is going on, if i’m writing

Left = (Lr + Li) * ds_excited

With the above code, do i “say” to Fenics that this expression is only evaluated on the ds_excited boundary?
When i’m trying to plot my solution, the following error occus:

fig1 = plt.figure()
fig11 = df.plot(p_r, title='part', cmap='jet')
fig1.colorbar(orientation="horizontal", mappable=fig11)
plt.show()

RuntimeWarning: invalid value encountered in reduce
  return umr_maximum(a, axis, None, out, keepdims, initial)
RuntimeWarning: invalid value encountered in less
  xa[xa < 0] = -1
RuntimeWarning: invalid value encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)

Apart from this my code runs perfectly fine, thank you very much!!

4

Hi, You don’t need the OnBoundary-class. You are also right that Left = (Lr + Li) * ds_excited is an expression only evaluated on the ds_excited boundary.

As you haven’t supplied your full code, I can’t tell why your plot doesn’t work. To debug this, I would first try to plot without a colorbar or save the function to pvd and visualize it with paraview.

1 Like
5

Plotting without a colorbar produced the same error.
I think i found my error, my solution variable only contains 0, i think it has to do with my mesh.
Thank your very much anyways, it really helped me out a lot. I just have to figure out how to do my mesh right.

6

Hi again,
theres one more question and answering this would help me out a lot.
When defining in gmsh

$PhysicalNames
4
2 1 "ExcitedBoundary"
2 2 "SoundProofBoundary"
2 3 "AdmittanceBoundary"
3 4 "DuctInside"
$EndPhysicalNames

for the interior domain for the model DuctInside, how do i use this information in my code?
Do I use

left_side = <expression> * df.dx

or

dx_volume = df.Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=4)
left_side = <expression> * dx_volume

Thanks very much for any help!!

7

Hi,
You should use the custom measures, as you sketch out.

1 Like
8

Thanks for the quick answer!
My code isn’t running, i always get NaN as a result, maybe you can have a quick view and tell me what i’m doing wrong because i really have no idea.
Gmsh-file:

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
4
2 1 "ExcitedBoundary"
2 2 "SoundProofBoundary"
2 3 "AdmittanceBoundary"
3 4 "DuctInside"
$EndPhysicalNames
$Nodes
9
1 1.7 0.02 -0.25
2 -1.7 0.02 -0.25
3 -1.7 -0.02 -0.25
4 1.7 -0.02 -0.25
5 1.7 0.02 0.25
6 -1.7 0.02 0.25
7 -1.7 -0.02 0.25
8 1.7 -0.02 0.25
9 0 0 0
$EndNodes
$Elements
24
1 2 2 1 1 2 7 3 
2 2 2 1 1 2 6 7
3 2 2 2 2 8 3 4
4 2 2 2 2 8 7 3
5 2 2 2 3 1 2 6
6 2 2 2 3 1 6 5
7 2 2 2 4 6 7 5
8 2 2 2 4 7 8 5
9 2 2 2 5 1 4 2
10 2 2 2 5 4 3 2
11 2 2 3 6 8 1 4
12 2 2 3 6 8 5 1
13 4 2 4 1 2 7 3 9
14 4 2 4 1 1 4 8 9
15 4 2 4 1 3 8 4 9
16 4 2 4 1 7 8 3 9
17 4 2 4 1 1 2 4 9
18 4 2 4 1 2 3 4 9
19 4 2 4 1 8 5 9 7
20 4 2 4 1 9 5 8 1
21 4 2 4 1 6 2 9 7
22 4 2 4 1 9 2 6 1
23 4 2 4 1 5 6 9 7
24 4 2 4 1 9 6 5 1
$EndElements

Code:

pathname = "./DuctMesh/"
msh = meshio.read(pathname + "DuctMesh.msh")

meshio.write(pathname + "DuctMesh_xdmf.xdmf", meshio.Mesh(points=msh.points,
                                                          cells={"triangle": msh.cells["triangle"]}
                                                          cell_data={"tetra": {"name_to_read": msh.cell_data["tetra"]["gmsh:physical"]}}))

meshio.write(pathname + "DuctMesh_xdmf_Boundary.xdmf",
             meshio.Mesh(points=msh.points,
                         cells={"triangle": msh.cells["triangle"]},
                         cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))

mesh = df.Mesh()

with df.XDMFFile(pathname + "DuctMesh_xdmf.xdmf") as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection("size_t", mesh, 2)
print("The MeshValueCollection: ", mvc)
with df.XDMFFile(pathname + "DuctMesh_xdmf_Boundary.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
print("The MeshFunctionSizet: ", mf)
# boundary = Bm.OnBoundary()  # define boundary
ds_excited = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=1)
ds_sound_proof = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_admittance = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)
dx_volume = df.Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=4)

# Create function space for complex variable
Vr = df.FiniteElement('Lagrange', mesh.ufl_cell(), 1)  # function space for real parts
Vi = df.FiniteElement('Lagrange', mesh.ufl_cell(), 1)  # function space for imag parts
Vc = df.FunctionSpace(mesh, Vr * Vi)  # composed function space for complex variable coefficients

# Define variational problem
(p_r, p_i) = df.TrialFunction(Vc)
(w_r, w_i) = df.TestFunction(Vc)

v_s = -df.Constant(1.0)  # - sign stems from direction of normal vector on excited boundary

# Define weak form - term by term
# stiffness matrix terms
t1r = df.inner(df.grad(p_r), df.grad(w_r)) - df.inner(df.grad(p_i), df.grad(w_i))
t1i = df.inner(df.grad(p_i), df.grad(w_r)) + df.inner(df.grad(p_r), df.grad(w_i))

# mass matrix terms
t3r = -(df.Constant(k) ** 2 * (df.inner(p_r, w_r) - df.inner(p_i, w_i)))
t3i = -(df.Constant(k) ** 2 * (df.inner(p_i, w_r) + df.inner(p_r, w_i)))

# Summation of matrix terms
ar = t1r + t3r
ai = t1i + t3i

# excitation term - best solution
Lr = omg*rho_f*v_s*w_i  # -omg*rho_f*v_s*w_r
Li = -omg*rho_f*v_s*w_r  # -omg*rho_f*v_s*w_i

a = (ar + ai) * dx_volume 
right= (Lr + Li) * ds_excited

p = df.Function(Vc)
df.solve(a == right, p)
(p_r, p_i) = p.split()

p_real = p_r.compute_vertex_values(mesh)
print(p_real)

Output: [nan nan nan nan nan nan nan nan  0.]

I would search by myself but there really is no proper documentation on how to do it or i’m not able to find something like that.
But thank you very much you already helped me out a lot!!

1 Like
9

Hi again,
can you maybe help me with my code i postet below? I really don’t know what to do, it just isn’t working.
Sorry to bother you that much.

10

Hi,

I’m having an issue converting a .msh file using the code snippets from above:

import meshio
import numpy as np
msh = meshio.read("graded_tri.msh")

line_cells = []
triangle_cells = []
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "line":
        if len(line_cells) == 0:
            line_cells = cell.data
        else:
            line_cells = np.vstack([line_cells, cell.data])

line_data = []
triangle_data = []
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        if len(line_data) == 0:
            line_data = msh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]])
    elif key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

triangle_mesh = meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells})
line_mesh =meshio.Mesh(points=msh.points,
                           cells=[("line", line_cells)],
                           cell_data={"name_to_read":[line_data]})
meshio.write("mesh.xdmf", triangle_mesh)

meshio.xdmf.write("mf.xdmf", line_mesh)

Gmsh version 4.6 was used to create and convert the following .geo file:
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {1, 1, 0, 1.0};
//+
Point(4) = {0, 1, 0, 1.0};
//+
Point(5) = {0.9, 0, 0, 1.0};
//+
Point(6) = {0.9, 0.9, 0, 1.0};
//+
Point(7) = {0.9, 1, 0, 1.0};
//+
Point(8) = {0, 0.9, 0, 1.0};
//+
Point(9) = {1, 0.9, 0, 1.0};
//+
Line(1) = {1, 5};
//+
Line(2) = {5, 2};
//+
Line(3) = {2, 9};
//+
Line(4) = {9, 3};
//+
Line(5) = {3, 7};
//+
Line(6) = {7, 4};
//+
Line(7) = {4, 8};
//+
Line(8) = {8, 1};
//+
Line(9) = {5, 6};
//+
Line(10) = {6, 8};
//+
Line(11) = {9, 6};
//+
Line(12) = {6, 7};
//+
Curve Loop(1) = {8, 1, 9, 10};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {9, -11, -3, -2};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {11, 12, -5, -4};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {7, -10, 12, 6};
//+
Plane Surface(4) = {4};
//+
Physical Curve(4) = {1};
//+
Physical Curve(5) = {2};
//+
Physical Curve(6) = {3};
//+
Physical Curve(7) = {4};
//+
Physical Curve(8) = {5};
//+
Physical Curve(9) = {6};
//+
Physical Curve(10) = {7};
//+
Physical Curve(11) = {8};

Using the version of meshio that was available from the repository about a week ago, I get the error message:

Traceback (most recent call last):
File “conversion.py”, line 31, in
meshio.write(“mesh.xdmf”, triangle_mesh)
File “…/env/lib/python3.6/site-packages/meshio/_helpers.py”, line 136, in write
if value.shape[1] != num_nodes_per_cell[key]:
AttributeError: ‘list’ object has no attribute ‘shape’

Please advise.

11

You have not added any physical surfaces to your geometry, thus meshio cannot find any.
Everything works if you add:

Physical Surface(1) = {1, 4, 3, 2};
12

You are a life saver dokken, thank you!

13

mesh
mesh1578×1078 20 KB

Hello again Dokken,

I have been revisiting this mesh and the conversion works beautifully. However, the actual xdmf files are lacking some information. As you can see in the mesh, I only have a single surface of 2d elements instead of four as expected.
It does not make a whole lot of sense to a gmesh noob.

Thanks,
eiriken

14

Please either attach a minimal example, or have a look at: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #104 by dokken which explain how to extract multiple regions. I have also made a tutorial on the subject: Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken

15

Certainly,
The following is the .geo file used to create the .msh file which was successfully converted using meshio convert:
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {1, 1, 0, 1.0};
//+
Point(4) = {0, 1, 0, 1.0};
//+
Point(5) = {0.9, 0, 0, 1.0};
//+
Point(6) = {0.9, 0.9, 0, 1.0};
//+
Point(7) = {0.9, 1, 0, 1.0};
//+
Point(8) = {0, 0.9, 0, 1.0};
//+
Point(9) = {1, 0.9, 0, 1.0};
//+
Line(1) = {1, 5};
//+
Line(2) = {5, 2};
//+
Line(3) = {2, 9};
//+
Line(4) = {9, 3};
//+
Line(5) = {3, 7};
//+
Line(6) = {7, 4};
//+
Line(7) = {4, 8};
//+
Line(8) = {8, 1};
//+
Line(9) = {5, 6};
//+
Line(10) = {6, 8};
//+
Line(11) = {9, 6};
//+
Line(12) = {6, 7};
//+
Curve Loop(1) = {8, 1, 9, 10};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {9, -11, -3, -2};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {11, 12, -5, -4};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {7, -10, 12, 6};
//+
Plane Surface(4) = {4};
//+
Physical Curve(4) = {1};
//+
Physical Curve(5) = {2};
//+
Physical Curve(6) = {3};
//+
Physical Curve(7) = {4};
//+
Physical Curve(8) = {5};
//+
Physical Curve(9) = {6};
//+
Physical Curve(10) = {7};
//+
Physical Curve(11) = {8};
//+
Physical Surface(1) = {1, 4, 3, 2};

The outputs are two .xdmf files, one with the solid portion shown in grey above and the mesh skeleton in colors.

It seems that only a portion of the physical surface is read to be 2d elements.

16

As I said in the previous post, follow my newest tutorial.
I copied your geo output and created an msh file (mesh2.msh), which I then converted with the following script from the tutorial:

import meshio
mesh_from_file = meshio.read("mesh2.msh")
import numpy

def create_mesh(mesh, cell_type, prune_z=False):
    cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
    cell_data = numpy.hstack([mesh.cell_data_dict["gmsh:physical"][key]
                         for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])
    # Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
    points= mesh.points
    if prune_z:
        points = points[:,:2]
    mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return mesh_new
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

which gives you the full triangular mesh

1 Like
17

Hello. I have a question. After creating mesh, I want to get information from this mesh and write it in .xdmf files. So, this is a code

import meshio
msh = meshio.read("/home/natalia/mesh/mesh_12.msh")
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "tetra":
        tetra_cells = cell.data

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write("mes_20.xdmf", tetra_mesh)

meshio.write("mf_20.xdmf", triangle_mesh)

I’m successfully get these files. After this I’m using Dolfin 2019.2.0.dev0 for creating .xml files which I can use in micromagnetic simulation, but Dolfin can’t recognise cell type. Reason: Unknown value “”. In line 27 of code in infile.read(mvc_tetra, “mes_20_name_to_read”).

from dolfin import *

# Create mesh from XDMF file
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

# Create a CellFunction to mark cell types
cell_markers = MeshFunction("size_t", mesh, 3)

# Iterate over cells and mark their types
for cell in cells(mesh):
    if CellType == "tetrahedron":
        cell_markers[cell] = 1  # tetrahedra are marked with value 1
    elif CellType == "triangle":
        cell_markers[cell] = 2  # triangles are marked with value 2

# Save cell information to XML file
with XDMFFile("cell_info.xml") as outfile:
    outfile.write(cell_markers)

# Create subdomain cell function for tetrahedra
mvc_tetra = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mes_20.xdmf") as infile:
    infile.read(mvc_tetra, "mes_20_name_to_read")

subdomain_marker_tetra = cpp.mesh.MeshFunctionSizet(mesh, mvc_tetra)

# Save subdomain information to XML file
with XDMFFile("subdomain_info_tetra.xml") as outfile:
    outfile.write(subdomain_marker_tetra)

# Create boundary subdomain facet function for triangles
mvc_triangle = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf_20.xdmf") as infile:
    infile.read(mvc_triangle, "mf_20_name_to_read")

boundary_marker_triangle = cpp.mesh.MeshFunctionSizet(mesh, mvc_triangle)

# Save boundary information to XML file
with XDMFFile("boundary_info_triangle.xml") as outfile:
    outfile.write(boundary_marker_triangle)

When I open xdmf files I can see, that they contain information about not cell but Topology. This is the whole data which file has.

<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XYZ"><DataItem DataType="Float" Dimensions="13195 3" Format="HDF" Precision="8">mes_20.h5:/data0</DataItem></Geometry><Topology TopologyType="Tetrahedron" NumberOfElements="63663" NodesPerElement="4"><DataItem DataType="Int" Dimensions="63663 4" Format="HDF" Precision="4">mes_20.h5:/data1</DataItem></Topology></Grid></Domain></Xdmf>

How Dolfin can recognise cell type in my case then?

18

There is topology information in the file.
However, you should not use the XML extension for xdmf files:

19

I use XML file for saving (outfile). Is it impossible? For taking information I take file XDMF:

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

# Create subdomain cell function for tetrahedra
mvc_tetra = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mes_20.xdmf") as infile:
    infile.read(mvc_tetra, "mes_20_name_to_read")

# Create boundary subdomain facet function for triangles
mvc_triangle = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf_20.xdmf") as infile:
    infile.read(mvc_triangle, "mf_20_name_to_read")
20

What I’m trying to tell you is that you cannot use the xml extension when writing an XDMFFile

should be

with XDMFFile("subdomain_info_tetra.xdmf") as outfile:
    outfile.write(subdomain_marker_tetra)

and
"mes_20_name_to_read" should just be "name_to_read"