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.