Import a mesh in FEniCS
Motivation and context
- When I first approached the software library FEniCS I was overwhelmed by the possibilities. This collection of software allows the user to define whatever Partial Differential Equation on a geometrical domain and automate all the process from the definition of finite elements to solution of the equation, it’s simply brilliant! Wanna solve a thermal diffusion problem? Just type the heat equation in integral form:
ET = ( rho*cspec*(T_trl-T_old)/dt*T_tst + \
+ Kcond*fe.dot(fe.grad(T_trl), fe.grad(T_tst)) )*fe.dx
- Wanna add convection through one of the boundaries? Just add that part to the equation:
ET = ( rho*cspec*(T_trl-T_old)/dt*T_tst +\
+ Kcond*fe.dot(fe.grad(T_trl), fe.grad(T_tst)) )*fe.dx +\
+ ( outer_hconv*(T_trl - outer_Temp)*T_tst )*ds_outer
- But here comes the catch, how to indicate to the software what part of the domain is
ds_outer? Seems like this is a limitation of the script like interface in FEniCS, using the provided functionalities (such as CompiledSubDomain), you can only define the geometry if you can type the equation for it.
- Well I refused to stop there, so I spent too many hours googling for a solution and I ended up collection enough pieces of code from the internet (mainly from Jørgen S. Dokken) to reach my goal: import mesh boundary definition from an external pre-processor to FEniCS.
Challenges I faced
- Before coming to FEniCS, the input mesh, coming from a pre-processor, has to be translated in a suitable format, I chose
xdmf and this can be achieved using the handy library meshio.
- In FEniCS a portion of the domain or of the boundary can be identified, or “marked”, with an integer identifier and called for the definition of a subset with specific material or a boundary condition.
- Before the “marking”, domain boundary and the materials subsets have to be created and it’s also a good idea to pre-set them to a null value:
# import the mesh in FEniCS
with fe.XDMFFile( meshfolder + meshfilename + ".xdmf") as xdmf:
mesh = fe.Mesh()
xdmf.read(mesh)
mesh.init()
# Define the whole boundary
boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
# initialize marker for the whole boundary
boundaries.set_all(0)
# define the whole domain
materials = fe.MeshFunction("size_t", mesh, mesh.topology().dim())
# initialize marker for the whole domain
materials.set_all(0)
- The geometrical subdomains can then be marked by looping over the cells and using the correspondence between the cell identification number and its property… which we still need to assign!
- The assignment has to be made in parallel to the translation operation from whatever format to
xdmf and stored in a separate object, a matrix is a good candidate, to be called when needed.
- Once all of this is done we get a nice table with the correspondence between the sets we defined in the pre-processor and the integers we can use in FEniCS for identification of geometrical features:
# import the mesh
mesh, boundaries, materials = importMesh.readMesh(meshfolder, meshfile)
# create the boundary integral measure
ds_inner = fe.Measure("ds", subdomain_id=1, subdomain_data=boundaries)
Useful lessons I learned
- This project didn’t involve much math or scientific thinking but it was critical to allow me to play with real parts coming from CAD into FEniCS!
- That’s it! I hope this information would turn out useful to anyone, I add the script I use for the conversion at the end of this page.
# -----------------------------------------------------------------70
# Import necessary packages
# -----------------------------------------------------------------70
import fenics as fe
import numpy as np
import meshio
import os
#%% ---------------------------------------------------------------70
# Auxiliary functions
# -----------------------------------------------------------------70
def fileprint(string, meshfolder, meshfilename):
print(string)
with open( meshfolder + meshfilename + ".md", "a") as file:
file.writelines( "\n" + string)
def getMeshInfo(meshfolder, meshfile, identifier = "cell", cell_geom = "triangle"):
"""Import the mesh from pre-processor and convert the mesh with Meshio
- use a specific keywork to identify cell sets in the mesh creation
and pass this keywork to the function as identifier
- specify whether to parse for triangle or tetra:
- 2D case --> cell_geom = "triangle"
- 3D case --> cell_geom = "tetra"
reference:
https://fenicsproject.discourse.group/t/definition-of-subdomain-by-element-sets/6440
"""
# Take the mesh file name without extension
meshfilename = meshfile.rsplit(".")[0]
# read the mesh whatever the format
pre_mesh = meshio.read(meshfolder + meshfile)
# import the cells sets
cell_sets = pre_mesh.cell_sets_dict
# import the cells
pre_cells = pre_mesh.get_cells_type(cell_geom)
# number of cells
n_cells = len(pre_cells)
# initialize a matrix to store sets marks for cell data (materials)
cell_data = np.zeros(n_cells)
# initialize a matrix to store sets marks for facet data (boundaries)
boundaries_data = np.zeros(n_cells)
# initialize marker counter
mrk = 0
# initialize marker dictionary
mrk_dict = {}
# loop through the dictionary items of cell_sets
# the key is the set name, the value is the sets cells (type and coordinates)
for marker, set in sorted(cell_sets.items()):
if "mesh" not in marker:
# Increse the marker for every set
mrk += 1
# Store set name and corresponding marker
mrk_dict[marker] = mrk
# loop through the types and values for a specific key
for type, entities in set.items():
# check for cell type
if type == cell_geom:
# Store in cell_data a marker with increasing int to be visualized in Paraview
boundaries_data[entities] = int(mrk)
# identify cell markers and store them in a list
if identifier in marker:
cell_data[entities] = int(mrk)
return pre_mesh, cell_data, boundaries_data, mrk, mrk_dict
def convertMesh(meshfolder, meshfile, identifier = "cell", cell_geom = "triangle"):
(pre_mesh, cell_data, boundaries_data, mrk, mrk_dict) = getMeshInfo(meshfolder, meshfile, identifier, cell_geom)
# import the cells
pre_cells = pre_mesh.get_cells_type(cell_geom)
# number of cells
n_cells = len(pre_cells)
# Take the mesh file name without extension
meshfilename = meshfile.rsplit(".")[0]
# Print the total number of cells
fileprint(f"Total number of cells in mesh = {n_cells}", meshfolder, meshfilename)
# Print the total number of markers
fileprint(f"Max marker index: {mrk:2d}", meshfolder, meshfilename)
# Print set name and corresponding marker
for key in mrk_dict.keys():
fileprint(f"Set: {key} -> marker {mrk_dict[key]}", meshfolder, meshfilename)
# Translate the mesh in xdmf format bringing the marker info
post_mesh = meshio.Mesh(points=pre_mesh.points,
cells={cell_geom: pre_cells},
cell_data={"name_to_read": [cell_data]} )
# Take the mesh file name without extension
meshfilename = meshfile.rsplit(".")[0]
# Adjust the filename and output the mesh file in xdmf format
meshio.write( meshfolder + meshfilename + ".xdmf", post_mesh)
def readMesh(meshfolder, meshfile, identifier = "cell", cell_geom = "triangle"):
# -----------------------------------------------------------------70
# Convert solid Element sets in surface boundary Facets Sets
# reference
# https://fenicsproject.discourse.group/t/how-to-define-boundary-conditions-on-an-irregular-geometry/2240/18
# -----------------------------------------------------------------70
# Take the mesh file name without extension
meshfilename = meshfile.rsplit(".")[0]
# check if mesh has already been converted in xdmf
xdmf_exists_flag = 0
for file in os.listdir(meshfolder):
if file == meshfilename + ".xdmf":
xdmf_exists_flag = 1
# convert the mesh if necessary
if xdmf_exists_flag == 0:
convertMesh(meshfolder, meshfile, cell_geom = cell_geom)
# import the mesh in FEniCS
with fe.XDMFFile( meshfolder + meshfilename + ".xdmf") as xdmf:
mesh = fe.Mesh()
xdmf.read(mesh)
mesh.init()
# Define the whole boundary
boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
# initialize marker for the whole boundary
boundaries.set_all(0)
# define the whole domain
materials = fe.MeshFunction("size_t", mesh, mesh.topology().dim())
# initialize marker for the whole domain
materials.set_all(0)
# initialize mesh to ask for facets
cell_to_facet = mesh.topology()(mesh.topology().dim(), mesh.topology().dim()-1)
# get cell_data from meshio
(pre_mesh, cell_data, boundaries_data, mrk, mrk_dict) = getMeshInfo(meshfolder, meshfile, identifier, cell_geom)
# loop through cells
for cell in fe.cells(mesh):
# mark the cell with info stored in cell_data
materials.set_value(cell.index(), int(cell_data[cell.index()]))
# loop through the facets of a cell
for face in cell_to_facet(cell.index()):
# extract Facet info with Dolfin internal method
facet = fe.Facet(mesh, face)
# mark the facet with the info stored in boundaries_data
if facet.exterior():
boundaries.set_value(facet.index(), int(boundaries_data[cell.index()]))
# https://www.karlin.mff.cuni.cz/~hron/fenics-tutorial/elasticity/doc.html
return mesh, boundaries, materials
#%% ---------------------------------------------------------------70
# Debugging check
# -----------------------------------------------------------------70
if __name__ == "__main__":
meshfolder = "./Meshes/"
meshfile = "WST.inp"
mesh, boundaries, materials = readMesh(meshfolder, meshfile, cell_geom = "triangle")
fe.XDMFFile( meshfolder + "check_boundaries.xdmf" ).write(boundaries)
fe.XDMFFile( meshfolder + "check_materials.xdmf" ).write(materials)
