Complex PDE on GMSH generated mesh

1

Hi everyone,
I’m currently trying to compute some acoustical scattering in a 3D channel with 1 subdomain and 3 different types of boundaries (2 Neumann and 1 Robin).
I have a costum generated mesh in gmsh using verison 4.3.0:
Gmsh-file:

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
4
2 1 "ExcitedBoundary"
2 2 "SoundProofBoundary"
2 3 "AdmittanceBoundary"
3 4 "Domain"
$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

Edit: For FEniCS i’m using verision 2019.1.0.
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)

with df.XDMFFile(pathname + "DuctMesh_xdmf_Boundary.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)


# define measure for excited boundary
ds_excited = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=1)
# define measure for sound proof boundary
ds_sound_proof = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
# define measure for admittance boundary
ds_admittance = df.Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)
# define measure for the inside domain
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 already tried the solutions found in this thread, but nothing of this works for me.
I just don’t understand why fenics only gives me NaN as an output.

Thanks in advance!

2

Hi,
It is simply because you havent defined a 3D mesh when writing it to XDMF.
It should be:

pathname = "./"
import meshio

msh = meshio.read(pathname + "juliandwain.msh")

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

In your code, cells where defined as triangle, thus you obtain a 2D surface mesh of the boundary of your channel.

1 Like
3

Hi, thanks for your answer!

Unfortunately, when I’m writing

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

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

as everything else stays the same, i get the following output

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:
[-inf  nan  nan  nan -inf -inf -inf  nan  nan]

I’m using the mesh from above and nothing else changed except “triangle” to “tetra”.

Thank you anyway for your answer, i really appreciate it!

4

Some notes, your dx_volume is not defined correctly. You should at least load the cell function (you are using a facet function in the measure).

mvc = df.MeshValueCollection("size_t", mesh, 3)
with df.XDMFFile(pathname + "juliandwain_xdmf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)

and
dx_volume = df.Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=4).
With your current code, you are not integrating over any volume. (As the volume is the full channel, you could just use dolfin.dx instead of defining a custom measure.

1 Like
5

Wow thank you very much, it works perfectly fine, you literally made my day!
I dind’t know i can use the dx command for gmsh generated meshes.

Just one note, if i’m trying to solve domains with two or more subdomains.

Do i write

with df.XDMFFile(pathname + "DuctMesh_xdmf_precise.xdmf") as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection("size_t", mesh)
with df.XDMFFile(pathname + "DuctMesh_xdmf_Boundary_precise.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc_2 = df.MeshValueCollection("size_t", mesh, 3)
with df.XDMFFile(pathname + "DuctMesh_xdmf_precise.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc_2)

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=cf, subdomain_id=4)

or

with df.XDMFFile(pathname + "DuctMesh_xdmf_precise.xdmf") as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection("size_t", mesh, 3)
with df.XDMFFile(pathname + "DuctMesh_xdmf_Boundary_precise.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
with df.XDMFFile(pathname + "DuctMesh_xdmf_precise.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)

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=cf, subdomain_id=4)

Thank your very much!!
You really helped me out a lot!

Edit 1:
In

mvc = df.MeshValueCollection("size_t", mesh, 3)

the last integer denotes the dimension, i guess, so if i’m trying to solve a 2D problem i would replace it with 2, wouldn’t I?

6

Hi, you need a separate MeshValueCollection for your facet data and your cell data. Thus, as described in the first snippet. The integer denotes the mesh entity dimension for the mesh value collection. Therefore, for a 2D mesh a cell function would have integer=2 and a facet function integer=1. For a 3D problem, a cell function has entity integer=3, facets have 2.

1 Like
7

Thank you very much again, I really appreciate it!
You literally solved all my problems I had the past weeks!