Subdomain for different materials based on physical volume

3

@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 = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data

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 = meshio.read(mshfile)

    ## 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 = cell.data
            else:
                tetra_cells = np.concatenate((tetra_cells,cell.data))
        elif cell.type == "triangle":
            if triangle_cells.all() == None:
                triangle_cells = cell.data
            else:
                triangle_cells = np.concatenate((triangle_cells,cell.data))

    ## put them together
    tetra_mesh = meshio.Mesh(points=msh.points, 
                             cells={"tetra": tetra_cells},
                             cell_data={"name_to_read":[tetra_data]})
    triangle_mesh =meshio.Mesh(points=msh.points,
                               cells=[("triangle", triangle_cells)],
                               cell_data={"name_to_read":[triangle_data]})
    
    ## 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:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "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: 
        infile.read(mesh)
        
    ## Import material info (physical volume)
    mvc = fn.MeshValueCollection("size_t", mesh, 2)
    with fn.XDMFFile("mesh.xdmf") as infile: 
        infile.read(mvc, "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: 
        infile.read(mvc, "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):
        super().__init__(**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
convertMesh('meshfilename.msh')
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