Hi, 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 https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html). 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.

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

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

Subdomain for different materials based on physical volume

Uncategorized

dokken June 16, 2020, 1:11am 2

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

Topic		Replies	Views
How to define different materials / importet 3D geometry (gmsh) Uncategorized	14	5519	November 4, 2020
How to use gmsh tags to assign material properties? mesh	7	1942	September 13, 2022
Assigning different physical groups to different color in gmsh and mark them as different material in fenics General	34	3636	June 6, 2023
Two materials and Mechanism General	11	923	June 27, 2023
How to import xdmf mesh generated by gmsh into fenics General	27	3145	January 20, 2023

Subdomain for different materials based on physical volume

Related topics