Mark regions from mesh of Gmsh

mesh
1

My problem is similar to Extracting mesh regions for gmsh mesh, except that I don’t solely focus on applying boundary conditions on surfaces.

My problem is that I have several 3D regions, each with different material properties (such as kappa=k0 in region 0, kappa=k1 in region 1, etc.)
When the mesh is created with mshr, there is already a working example in the tutorial and there’s no problem. However when the mesh is created with Gmsh and imported in FEniCS, I do not know how to mark the regions.

Here are some parts of my code:

k_0 = 18
k_1 = 34
k_2 = 10

# Define the values of kappa in each subdomain
class K(UserExpression):
    def __init__(self, materials, k_0, k_1, k_2, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.k_0 = k_0
        self.k_1 = k_1
        self.k_2 = k_2

def eval_cell(self, values, x, cell):
    if self.materials[cell.index] == 0:
        values[0] = self.k_0
    elif self.materials[cell.index] == 1:
        values[0] = self.k_1
    else:
        values[0] = self.k_2

def value_shape(self):
    return ()
kappa = K(materials, k_0, k_1, k_2)

And in the weak form there is “kappa”, so kappa takes its corresponding value according to the region the equation is being solved.

But for this snippet to work, 3 regions need to be marked as 0, 1 and 2. And this is what I struggle with.

I already have a “cf” and a “mf” defined as in Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio. I just do not know how to take it from there…

2

Just import the mesh and then create different domains. I had a mesh in .XML format and then I imported that into Fenics and created different domains like this.

import mesh

mesh = Mesh(‘halfcylinder.xml’)
gdim = mesh.geometry().dim()

mark the boundary

boundary_parts = MeshFunction(‘size_t’, mesh, mesh.topology().dim()-1)
y0 = AutoSubDomain(lambda x: x[0] < -0.5 and near(x[1], 0) )
y0.mark(boundary_parts, 1)

3

I would rather avoid dealing with numerical floats like the plague as it is source of errors and would require a manual adjustment anytime I modify the mesh, which is horrible and sounds like a badly thought hack.

I already have 3 distinct regions, and defining volume measures with them is easy:
dx_custom0 = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=1)

I wish I could do something like

cf(id=1).mark(0)
But obviously this doesn’t work.

This would be insanely better than manually set the regions, which may not even correspond to how the mesh is defined…

4

I guess your are looking for something similar to: Dolfinx discontinous expression

As long as you can access the underlying PETScVector in the dolfin.Function, you can follow this thread