Extracting mesh regions for gmsh mesh

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!

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

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!!

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

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.

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!!

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

1 Like

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!!