Working with different subdomains: from meshio to dolfin

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.

Please mark the Physical Surface and Curve with numbers, not strings. As you can observe by inspecting the boundary.xdmf file the values on each boundary is not 1 and 2, but 2 and 3, and does not consist of the whole circle.
This is because you are converting all the data written in with meshio.
See for instance: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #104 by dokken
or Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken
for how to convert a meshio mesh and corresponding boundary data in a safe way.

1 Like

Thank you very much @dokken, with your advice I was able to get the correct results. Just for clarity, if anyone is interested, I’m posting the .geo and the .py corrected versions below. I have two last questions, one is about this statement:

What does it mean? What is the command that has changed from the previous to the final version?

The other question regards the "name_to_read" string that appears, for example, within meshio.Mesh() . Is it an arbitrary string or do we need to use this specific string?

Thank you very much again.


As promised, .geo file:

//+
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 Curve(1) = {1};
//+
Physical Curve(2) = {2};
//+
Physical Surface(3) = {1};

and .py file:

import meshio
import numpy as np

msh = meshio.read("donut-2d.msh")

cells = np.vstack(np.array([cells.data for cells in msh.cells
                            if cells.type=="triangle"]))

triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", cells)])

facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                  if cells.type=="line"]))

facet_data = msh.cell_data_dict["gmsh:physical"]["line"]

facet_mesh = meshio.Mesh(points=msh.points,
                         cells=[("line", facet_cells)],
                         cell_data={"name_to_read": [facet_data]})

meshio.write("boundaries.xdmf", facet_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("boundaries.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)

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(1)))
print("epsilon length: ", assemble(Constant(1)*ds(2)))

I’ve not tested the old code for a while now, so I would suggest you sit down and compare the inputs you are giving the meshio meshes.

You can use any string, as long as you are consistent with using the same string when reading in the mesh data.

1 Like

Hi,

I’d like to have an advice for gmsh. I think that it is not necessary to open another topic, so I’m writing here.

I’m trying to create a spherical shell in gmsh, with a mesh that is denser near the inner boundary of the shell.

For this purpose, I wrote the following .geo file, and then I executed the 3D command in the gmsh gui for generating the 3D mesh, that is then saved in a .msh file:

SetFactory("OpenCASCADE");
Sphere(1) = {0, 0, 0, 0.1, -Pi/2, Pi/2, 2*Pi};
//+
Sphere(2) = {0, 0, 0, 0.12, -Pi/2, Pi/2, 2*Pi};
//+
Sphere(3) = {0, 0, 0, 0.4, -Pi/2, Pi/2, 2*Pi};
//+
BooleanDifference(4) = { Volume{3};  Delete;}{ Volume{2};  };  
//+
BooleanDifference(5) = { Volume{2};  Delete;}{ Volume{1}; Delete; };
//+
BooleanFragments{ Volume{5}; Volume{4}; Delete; }{ }
//+
Characteristic Length{ PointsOf{ Volume{4}; } } = 0.04;
//+
Characteristic Length{ PointsOf{ Volume{5}; } } = 0.005;
//+
Physical Surface(1) = {4};
//+
Physical Surface(2) = {3};
//+
Physical Volume(3) = {4, 5}; 

I used the second sphere to create a dummy volume for imposing a denser mesh step. The generated mesh seems visually good to me:

In order to see if everything works correctly in FEniCS, I first check the volume and surface integral of the spherical shell, of the inner boundary and of the outer boundary respectively, with the following .py script:

from dolfin import *
import meshio
import numpy as np

def msh_to_mf(dirname):

    msh = meshio.read(dirname + "spherical-shell-3d.msh")

    cells = np.vstack(np.array([cells.data for cells in msh.cells
                                if cells.type=="tetra"]))

    tetra_mesh = meshio.Mesh(points=msh.points,
                                cells=[("tetra", cells)])

    facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                      if cells.type=="triangle"]))

    facet_data = msh.cell_data_dict["gmsh:physical"]["triangle"]

    facet_mesh = meshio.Mesh(points=msh.points,
                             cells=[("triangle", facet_cells)],
                             cell_data={"name_to_read": [facet_data]})

    meshio.write(dirname + "boundaries.xdmf", facet_mesh)
    meshio.write(dirname + "volume.xdmf", tetra_mesh)

    parameters["ghost_mode"] = "shared_vertex"
    mesh = Mesh()

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

    mvc_2d = MeshValueCollection("size_t", mesh, 2)

    with XDMFFile(dirname + "boundaries.xdmf") as infile:
        infile.read(mvc_2d, "name_to_read")

    mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)

    mvc_3d = MeshValueCollection("size_t", mesh, 3)

    mf_3d = cpp.mesh.MeshFunctionSizet(mesh, mvc_3d)

    return mesh, mf_2d, mf_3d

mesh, mf_2d, mf_3d = msh_to_mf("./")

ds = Measure("ds", domain=mesh, subdomain_data=mf_2d)
dx = Measure("dx", domain=mesh, subdomain_data=mf_3d)

print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*ds))
print(assemble(Constant(1)*ds(1)))
print(assemble(Constant(1)*ds(2)))
print(assemble(Constant(1)*ds(3)))

The integral seems to be also correct.

I’d like to know if there is a simpler and "fenics-friendly" way to write the .geo file, since, even if it seems to be correct, I’m afraid that when I’m going to run the simulation there will be some errors arising due to the cumbersomeness of the .geo file.