Subdomain for different materials based on physical volume


I hope there is an existing answer:
I have a complicated geometry with multiple materials combined, and I wish to use constant values that are unique to each material (e.g. Kappa in

These materials are defined as Physical Volume in the mesh, so I wish to extract this data and apply that information to define these constants (k=1 for Physical Volume 1, k=1.25 for Physical Volume 2, etc.) without definind the boundary with mathematical expression.

I see many posts asking similar questions, but I wasn’t able to find the answer that works for my situation.
Can someone guide me to a discussion thread talking about this or provide a minimum working example to do this? Thank you.

1 Like

See How to define different materials / importet 3D geometry (gmsh)

1 Like

@dokken Thank you for the suggestion. I have solved my issue, but for other people who might be looking for a tip (that are not trivial to beginners), I will describe my issue and the solution below:

After hours of investigation, I have noticed there were three main issues I’ve had:

1: xdmf output

My model is consisted of two objects, and each of them are stored as different elements in msh.cells array. This means there are two separated cell data sections whose type is ‘tetra’ (and same goes to ‘triangle’), and simply doing

for cell in msh.cells:
    if cell.type == "tetra": 
        tetra_cells =
    elif cell.type == "triangle":
        triangle_cells =

will fail to include all of such cell data because this code overwrites every time it encounters a new cell data array. To fix this, I changed the code to cocatenate those data. Below is a minimum working example of my mesh convert function. I am sure there is more elegant way of doing this, but at least this code works.

def convertMesh(mshfile):
    ## original mesh file
    msh =

    ## physical surface & volume data
    for key in msh.cell_data_dict["gmsh:physical"].keys():
        if key == "triangle":
            triangle_data = msh.cell_data_dict["gmsh:physical"][key]
        elif key == "tetra":
            tetra_data = msh.cell_data_dict["gmsh:physical"][key]

    ## cell data
    tetra_cells = np.array([None])
    triangle_cells = np.array([None])
    for cell in msh.cells:
        if cell.type == "tetra":
            if tetra_cells.all() == None:
                tetra_cells =
                tetra_cells = np.concatenate((tetra_cells,
        elif cell.type == "triangle":
            if triangle_cells.all() == None:
                triangle_cells =
                triangle_cells = np.concatenate((triangle_cells,

    ## put them together
    tetra_mesh = meshio.Mesh(points=msh.points, 
                             cells={"tetra": tetra_cells},
    triangle_mesh =meshio.Mesh(points=msh.points,
                               cells=[("triangle", triangle_cells)],
    ## output
    meshio.write("mesh.xdmf", tetra_mesh)
    meshio.write("mf.xdmf", triangle_mesh)

2. Physical volume data import

As a beginner, it took me a while until I realize I was not importing the physical volume data in a code that is copied-and-pasted from the long thread about meshio. The often-quoted example

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

is only looking for the "name_to_read" in mf.xdmf, which returns (according to the convertMesh() above) triangle_data. This is, of course, only returning the physical surface data rather than the physical volume data which I want.
Although I am not 100% sure about the treatment of mvc, the function below is what I am using to import the mesh (and it works for me):

def importMesh():
    ## Import mesh
    mesh = fn.Mesh()
    with fn.XDMFFile("mesh.xdmf") as infile:
    ## Import material info (physical volume)
    mvc = fn.MeshValueCollection("size_t", mesh, 2)
    with fn.XDMFFile("mesh.xdmf") as infile:, "name_to_read")
    materials = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc)

    ## Import boundary info (physical surface)
    # overwriting mvc because it's not really being used anywhere else
    mvc = fn.MeshValueCollection("size_t", mesh, 2) 
    with fn.XDMFFile("mf.xdmf") as infile:, "name_to_read")
    boundaries = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    return mesh,materials,boundaries

3. Defining a constant (Coefficient) based on materials

This is what I wanted to do. The goal was to assign the density to each material (\rho: ‘rho’) that are marked individually with physical volume id. This density information rho is used when evaluating the PDE. To do this, the following class is used:

class Rho(fn.UserExpression):
    def __init__(self,materials,volume_list,rho_list,**kwargs):
        self.materials = materials
        self.volume_list = volume_list
        self.rho_list = rho_list
    def value_shape(self):
        return ()
    def eval_cell(self,values,x,cell):
        values[0] = 0
        for i in range(len(self.volume_list)):
            if self.materials[cell.index] == self.volume_list[i]:
                values[0] = self.rho_list[i]   

working example:

Using above functions & a class, my code looks like this:

import fenics as fn
import numpy as np
import meshio

# mesh prep
mesh,materials,boundaries = importMesh()

# function space
dx = fn.Measure('dx', domain=mesh, subdomain_data=mvc)
V = fn.FunctionSpace(mesh, 'CG', 1)

# boundary condition 
# (applying a pre-defined boundary value outer_value to physical surface #5)
outer_boundary = fn.DirichletBC(V, outer_value, boundaries, 5)
bcs = outer_boundary

# material assignment
# (applying three different values to physical volume #3,4, and 6)
volume_list = [3,4,6]                          # list of physical volume ID's
rho_list    = [rho_vacuum,rho_inner,rho_outer] # list of pre-defined values
rho = Rho(materials,volume_list,rho_list,degree=0)

# solve (just an example: not really relevant)
# solving a PDE for a scalar field phi (nonlinear Poisson)
# pre-defined constants: initial_guess, Lambda, M, scaling
phi = fn.project(initial_guess,V)
v = fn.TestFunction(V)
f = Lambda**5/phi**2
g = fn.Constant(0)
a = (fn.inner(fn.grad(phi)*scaling**2,fn.grad(v))+ rho*v/M)/scaling*fn.dx
L = f*v/scaling*fn.dx + g*v*fn.ds
F = a-L
fn.solve(F==0, phi, bcs)

I hope this helps someone who is stuck at the same place as I was.

1 Like

hello, Sterling! I am also tryng to define different materials. and your solution helps me a lot. I really appreciate it !