# Mark regions from mesh of Gmsh

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…

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)

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…

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