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 = self.k_0 elif self.materials[cell.index] == 1: values = self.k_1 else: values = 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…