@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.