Hi all,
I’m trying to understand the correct workflow for meshing in fenics, starting from a .msh file and arriving to a working Measure in fenics. In particular, I’m working on this example: my geometry consists of a disk with a hole at the center (poorly speaking, a 2D donut). I created this geometry in gmsh writing the following .geo script:
//+
SetFactory("OpenCASCADE");
Circle(1) = {0.5, 0.5, 0.5, 0.1, 0, 2*Pi};
//+
Circle(2) = {0.5, 0.5, 0.5, 1, 0, 2*Pi};
//+
Curve Loop(1) = {2};
//+
Curve Loop(2) = {1};
//+
Plane Surface(1) = {1, 2};
//+
Physical Surface("omega") = {1};
//+
Physical Curve("gamma") = {1};
//+
Physical Curve("epsilon") = {2};
I labeled the surface of the geometry with omega, while the inner and outer circles with gamma and epsilon, respectively.
Then, I obtained the .msh file by using the 2D and the Refine by splitting commands (refined two times) in the gmsh gui.
In order to read the mesh and import it in FEniCS, I’ve written the following Python script:
import meshio
msh = meshio.read("donut-2d.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "line":
line_data = msh.cell_data_dict["gmsh:physical"][key]
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
if cell.type == "line":
line_cells = cell.data
if cell.type == "triangle":
triangle_cells = cell.data
line_mesh = meshio.Mesh(points=msh.points,
cells=[("line", line_cells)],
cell_data={"name_to_read":[line_data]})
triangle_mesh = meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read":[triangle_data]})
meshio.write("boundary.xdmf", line_mesh)
meshio.write("surface.xdmf", triangle_mesh)
from dolfin import *
mesh = Mesh()
with XDMFFile("surface.xdmf") as infile:
infile.read(mesh)
mvc_1d = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("boundary.xdmf") as infile:
infile.read(mvc_1d, "name_to_read")
mf_1d = cpp.mesh.MeshFunctionSizet(mesh, mvc_1d)
mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("surface.xdmf") as infile:
infile.read(mvc_2d, "name_to_read")
mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
ds = Measure("ds", domain=mesh, subdomain_data=mf_1d)
dx = Measure("dx", domain=mesh, subdomain_data=mf_2d)
print("")
print("omega area: ", assemble(Constant(1)*dx))
print("gamma + epsilon length: ", assemble(Constant(1)*ds))
print("gamma length: ", assemble(Constant(1)*ds(2)))
print("epsilon length: ", assemble(Constant(1)*ds(3)))
The last lines give me the following outputs:
omega area: 3.1097676855865237
gamma + epsilon length: 6.910189733699893
gamma length: 1.0926747887983828
epsilon length: 5.190205246792324
Reading other posts, I interpreted the quantities within the print function as the integrals over the different measures (is it right?). As you can see, the integral over dx and ds are similar to the theoretical values, while this is not true for ds(2) and ds(3).
So I think I’m quite out of the way, and I’m not working correctly. What I want to do in dolfin is to integrate some terms over gamma and some other terms over epsilon in my weak formulation. Additionally, I could want to use one of these curves for boundary conditions.
So my question is: what is the correct way of passing and/or calling two different (one-dimensional in this case) subdomains from meshio to dolfin?
I’m aware that very probably I’ve misunderstood the meaning of some variables, and I have to admit that I’ve made a sort of “leap of faith” writing this script, since I do not get what kind of operations are done with meshio in the code above, so I’ve just put together what I learned from the posts in this group, but probably this was not enough.
Thank you in advance for your precious help.
