Mesh input/output

Introduction

This tutorial will demonstrate how you can import geometric computational meshes to use for your finite element simulations. For this demonstration we will be loading in a 3D simplex tetrahedral mesh of a section of a cardiomyocyte geometry.

Learning objectives

• Learn how to read in mesh files and use them to define the OpenCMISS-Iron objects that are related to mesh nodes and elements for a finite element problem.
• Know how meshes can be output into different formats.

Running this tutorial

Provided that you have already installed the OpenCMISS library, the code excerpts in this section can be run interactively by entering it directly into the Python interpreter, which can be started by running python or ipython in a terminal. The code shown here can also be downloaded as a jupyter notebook in the following link.

Dependencies

The meshio python module is required for running this tutorial. If you are running this tutorial from an OpenCMISS-Iron docker container, then this module is already installed.

Loading in a 3D tetrahedral mesh

In this tutorial we will walk through the steps to load in a 3D finite element linear simplex mesh of a section of a cardiomyocyte. The filenames are input_mesh.node and input_mesh.ele. The file extensions indicate that the .node file contains nodal information and the .ele file contains element information. The mesh files contain header information in the first line and then list the mesh node and element properties, respectively .The mesh we are using in this tutorial was originally generated using a software package for creating tetrahedral meshes, called tetgen but we have stripped away some of the details to maintain the focus of this tutorial on using the OpenCMISS-Iron libraries. ## Key concepts when loading in a mesh geometry

When you start considering importing meshes into an OpenCMISS-Iron finite element simulation keep the following points in mind:

1. Remember that OpenCMISS-Iron is just a library of functions to perform finite element simulations. Open-CMISS Iron does not have a pre-defined file format that it recognises. OpenCMISS-Iron is therefore mesh-file-format-agnostic. This is to say that you can import any type of file you have into the python program for your finite element simulation.
2. The key step to importing meshes into your program is to write some code that enables you to import the mesh nodes and elements into aarray. This array is then used to define the Nodes and Mesh Elements objects that define a finite element mesh. The mesh format inputted is not important so long as you can use python to read that file and convert the data to a numpy array of element matrix (connectivity matrix i.e. elements and their associated nodes) and a nodal matrix (nodes and their associated x-y coordinates)
3. When importing meshes you need to know what type of interpolation scheme is used and then define the corresponding basis functions within the finite element program that you are using for the finite element simulation.
4. Be very careful and aware of the node numbering sequence that you use to define each finite element. OpenCMISS-Iron uses the Zienkiewicz node numbering system for each element:

5. Output of meshes is the same as input - you just need to extract the node and element variables from the OpenCMISS-Iron objects and write them out in the format of your choice.
6. At the end of this tutorial we provide links to a mesh import and export functions that we have written previously, which you are welcome to use and customise for your own purposes.

Types of Basis Functions Compatible with OpenCMISS-Iron

Currently, the types of basis functions that have been defined in OpenCMISS-Iron include:

Linear Lagrange

Quadratic Lagrange

Cubic-Hermite

Linear Simplex

Quadratic Simplex

Loading the OpenCMISS-Iron library

In order to use OpenCMISS we have to first import the opencmiss.iron module from the opencmiss package.

# Intialise OpenCMISS-Iron.
from opencmiss.iron import iron

Assuming OpenCMISS has been correctly built with the Python bindings by following the instructions in the programmer documentation, we can now access all the OpenCMISS functions, classes and constants under the iron namespace.

The next section describes how we can interact with the OpenCMISS-Iron library through an object-oriented API.

Step by step guide

1. Creating a coordinate system

First we construct a coordinate system that will be used to describe the geometry in our problem. The 3D geometry will exist in a 3D space, so we need a 3D coordinate system.

# Create coordinate system.
coordinate_system_user_number = 1
coordinate_system = iron.CoordinateSystem()
coordinate_system.CreateStart(coordinate_system_user_number)
coordinate_system.DimensionSet(3)
coordinate_system.CreateFinish()

2. Read in the node file

OpenCMISS-Iron can handle any user-generated meshes provided the files are processed and entered into OpenCMISS-Iron objects related to the mesh. This process requires that the user provide files detailing the nodes and elements. We will start by reading the details of the nodes into our OpenCMISS-Iron based python program.

# Reading node file
node_file = open('input_mesh.node', 'r')

# Reading the mesh details from the first line of the node file
num_nodes, num_coords, num_attributes, boundary_markers = node_file.readline().split()
num_nodes = int(num_nodes)
num_coords=int(num_coords)
num_attributes=int(num_attributes)
boundary_markers=int(boundary_markers)
# Creating variables to store node number & boundary marker
NodeNums = [[0 for m in range(2)] for n in range(num_nodes)]
# Creating array to store x and y coordinates
NodeCoords = [[0 for m in range(num_coords)] for n in range(num_nodes)]

# Reading details from node file
for i in range(num_nodes):
    NodeNums[i][0], NodeCoords[i][0], NodeCoords[i][1], NodeCoords[i][2] = node_file.readline().split()  #node number, x, y, z, boundary marker
    # Converting from string to appropriate type
    NodeNums[i][0] = int(NodeNums[i][0])
    NodeCoords[i][0] = float(NodeCoords[i][0])
    NodeCoords[i][1] = float(NodeCoords[i][1])
    NodeCoords[i][2] = float(NodeCoords[i][2])

# Closing the file
node_file.close()

3. Read in the elements

OpenCMISS-Iron uses Zeinkiewicz winding to map out the nodes of each element. This means that we need to manually map each element’s nodes from the order in which they appear in the file to that required by Iron. In this example the order of the nodes in the file matches the OpenCMISS-Iron winding.

#Inputing element file
elem_file = open('input_mesh.ele', 'r')

#Reading the values of the first line
num_elements, nodes_per_ele, ele_attributes = elem_file.readline().split()

# Converting the numbers to integers
num_elements = int(num_elements)
nodes_per_ele = int(nodes_per_ele)
ele_attributes = int(ele_attributes)

# Creating variable to store the element map
ElemMap = [[0 for x in range(nodes_per_ele+ele_attributes)] for y in range(num_elements)]

Elemindex = [[0 for m in range(1)] for n in range(num_elements)]
# Reading element data from elemfile
for i in range(num_elements):
    # Performing the mapping
    Elemindex[i][0], ElemMap[i][0], ElemMap[i][1], ElemMap[i][2], ElemMap[i][3] = elem_file.readline().split()


# Closing the file
elem_file.close()

4. Creating basis functions

The finite element description of our fields requires a basis function to interpolate field values over elements, so we create a 3D basis with simplex interpolation in all \(\xi\) directions.

numberOfXi = 3
basis_user_number = 1
basis = iron.Basis()
basis.CreateStart(basis_user_number)
basis.type = iron.BasisTypes.SIMPLEX
basis.numberOfXi = numberOfXi
basis.interpolationXi = [iron.BasisInterpolationSpecifications.LINEAR_SIMPLEX] * numberOfXi
basis.quadratureOrder = 2
basis.CreateFinish()

5. Creating a region

Next we create a region that our fields will be defined on and tell it to use the 3D coordinate system we created previously. The CreateStart method for a region requires another region as a parameter. We use the world region that is created by default so that our region is a subregion of the world region.

# Create region
region_user_number = 1
region = iron.Region()
region.CreateStart(region_user_number, iron.WorldRegion)
region.CoordinateSystemSet(coordinate_system)
region.LabelSet("Region")
region.CreateFinish()

6. Define mesh objects

Finally, we need to intialise the nodes, the elements and the mesh with the mesh parameters that we have imported into the python program.

First, lets initialise an iron.Nodes() object which will store the details of the nodes of the mesh.

# Initialise Nodes
nodes = iron.Nodes()
nodes.CreateStart(region, num_nodes)
nodes.CreateFinish()

Then, initialise an iron.Mesh() object to store the details of the mesh itself.

# Initialise Mesh
mesh = iron.Mesh()
mesh_user_number=1
mesh.CreateStart(mesh_user_number, region, num_coords)
mesh.NumberOfElementsSet(num_elements)
mesh.NumberOfComponentsSet(1)

We now create a meshelements object which will define the node numbering for each element. Tetgen generates a .ele file that contains a sequence of nodes to define each element. As mentioned at the start of this tutorial OpenCMISS-Iron uses the Zienkiewicz winding but tetgen uses a different winding scheme for the node numbering. Therefore we have to ensure that we alter the node numbering when defining the mesh elements below to reflect the OpenCMISS-Iron winding.

# Initialise Elements
meshElements = iron.MeshElements()
meshElements.CreateStart(mesh, 1, basis)

for i in range(num_elements):
    element=Elemindex[i][0]
    nodesList = list(
      map(int,[ElemMap[i][0], ElemMap[i][1], ElemMap[i][2], ElemMap[i][3]]))
    meshElements.NodesSet(int(element), nodesList)

meshElements.CreateFinish()

We conclude this step by finalizing the mesh object:

# Finalise Mesh
mesh.CreateFinish()

7. Decomposing the mesh

Once the mesh has been created we can decompose it into a number of domains in order to allow for parallelism. We choose the options to let OpenCMISS calculate the best way to break up the mesh. We also set the number of domains to be equal to the number of computational nodes this example is running on. Note that if MPI infrastructure is not used, only single domain will be created. Look for our parallelisation example for an illustration of how to execute simulations using parallel processing techniques.

# Perform mesh decomposition.
decomposition_user_number = 1
decomposition = iron.Decomposition()
decomposition.CreateStart(decomposition_user_number, mesh)
decomposition.CreateFinish()

8. Creating a geometric field

Now that the mesh has been decomposed we are in a position to create fields. The first field we need to create is the geometric field. Here we create a field and partition the field to different computational nodes using the mesh decomposition that we have just created. Once we have finished creating the field we can change the field DOFs to define the geometry using the nodal coordinates we have imported.

# Create a field for the geometry.
geometric_field_user_number = 1
geometric_field = iron.Field()
geometric_field.CreateStart(geometric_field_user_number, region) #notice that the geometric field is associated with region in this function call.
geometric_field.LabelSet('Geometry')
geometric_field.MeshDecompositionSet(decomposition)
geometric_field.CreateFinish()

We have created a geometric field but all the field component values are currently set to zero by default. We can define the geometry for the custom mesh as follows. The key step is the call to the geometric field function ParameterSetUpdateNodeDP. You can find details about this method within the OpenCMISS-Iron python API. Since OpenCMISS-Iron has been coded with MPI parallel programming in mind we access the node variables within the decomposition object that has been defined so that we correctly update node variables in the right computational node. Of course, if you are only using 1 computational node then this code will revert to updating all the nodes in the mesh on that single computational node.

# Set geometric field values from customized mesh.
computationalNodeNumber = iron.ComputationalNodeNumberGet()
for node_idx in range(num_nodes):
    node_id = NodeNums[node_idx][0]
    node_x = NodeCoords[node_idx][0]
    node_y = NodeCoords[node_idx][1]
    node_z = NodeCoords[node_idx][2]
    geometric_field.ParameterSetUpdateNodeDP(
      iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES,
      1, 1, node_id, 1, node_x)
    geometric_field.ParameterSetUpdateNodeDP(
      iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES,
      1, 1, node_id, 2, node_y)
    geometric_field.ParameterSetUpdateNodeDP(
      iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES,
      1, 1, node_id, 3, node_z)

Once the geometric field variables have been updated with the input mesh node coordinates it is important to ensure that all computational nodes are notified with the new information. So the following function calls are really important.

# Update the geometric field.
geometric_field.ParameterSetUpdateStart(
  iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES)
geometric_field.ParameterSetUpdateFinish(
  iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES)

9. Exporting meshes

Now that the mesh has been imported into the OpenCMISS-Iron mesh and related objects, we can now use them for solving equations on the mesh topology. In what follows we show examples of how information in mesh objects can be extracted to generate mesh outputs. Typically, we don’t want to only visualise the mesh node coordinates and elements. We usually execute a biophysics simulation over the input mesh domain and it is the resulting solution over the mesh that we want to visualise. The solution is typically stored in a field that is created when setting up the simulation. We have not set up a simulation in this example. Therefore, for illustration purposes lets create an additional field that is associated with the mesh. We will then set out to visualise this “output” over the mesh in t he examples that follow.

### Setting up the part of the example about outputting mesh info.
## create new field
solution_field_user_number=2
solution_field=iron.Field()
#associate the field with the region and mesh decomposition we have set up in this example.
solution_field.CreateStart(solution_field_user_number,region)
solution_field.LabelSet('Solution')
solution_field.NumberOfComponentsSet(iron.FieldVariableTypes.U,1)
solution_field.MeshDecompositionSet(decomposition)
solution_field.CreateFinish()
#initialise the values in the field to some arbitrary value for this example.
solution_field.ComponentValuesInitialiseDP(
  iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,1000.0)

Example 1: Output as list of nodes and elements

For a given mesh, one can access the mesh nodes as follows:

# meshNodes = iron.MeshNodes()
# mesh.NodesGet(1,meshNodes)
# num_nodes= meshNodes.NumberOfNodesGet()

Since we already have nodes defined, lets just use that in this example. First, lets access the number of coordinates and nodes in the mesh:

#since we already have Nodes defined, lets just use thatin this example
num_nodes = nodes.NumberOfNodesGet()
num_coords = coordinate_system.DimensionGet()

Lets create a nodes list that contains the node numbers, the coordinates of each node and the solution value at each node. We use the iron function Field.ParameterSetGetNodeDP() to get the values from the geometric and solution fields.

nodes_list = [[0 for coord in range(0,num_coords+2)] for node in range(0,num_nodes)]


#create list object containing nodes and another one for elements
for i in range(0,num_nodes):
    node=i+1
    node_coordx = geometric_field.ParameterSetGetNodeDP(
      iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,1,node,1)
    node_coordy = geometric_field.ParameterSetGetNodeDP(
      iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,1,node,2)
    node_coordz = geometric_field.ParameterSetGetNodeDP(
      iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,1,node,3)
    node_solution_value = solution_field.ParameterSetGetNodeDP(
      iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,1,node,1)
    nodes_list[i] = [node,node_coordx,node_coordy,node_coordz,node_solution_value]

Write out a file containing the list we have just created.

node_file = open('output_mesh.node', 'w')
node_file.writelines([str(line) + "\n" for line in nodes_list])
# Closing the file
node_file.close()

Now lets create a list of elements and write them out similarly:

num_elem = mesh.NumberOfElementsGet()
elements_list = [[0 for i in range(5)] for elem in range(num_elem)]
for j in range(0,num_elem):
    elemidx=j+1
    elemnodes = meshElements.NodesGet(elemidx,4)
    elements_list[j] = [
      elemidx,
      elemnodes[0],
      elemnodes[1],
      elemnodes[2],
      elemnodes[3]]

elem_file=open('output_mesh.ele','w')
elem_file.writelines([str(line) + "\n" for line in elements_list])
elem_file.close()

Example 2: Output as numpy arrays - useful for importing into ParaView

#Example 2: Output as arrays
import numpy as np

nodesx = np.array(nodes_list[:][1])
np.save('output_mesh.x.npy',nodesx)
nodesy = np.array(nodes_list[:][2])
np.save('output_mesh.y.npy',nodesy)
nodesz= np.array(nodes_list[:][3])
np.save('output_mesh.z.npy',nodesz)
nodesolution=np.array(nodes_list[:][4])
np.save('output_mesh.sol.npy',nodesolution)

# Export results in Exfile format.
# fields = iron.Fields()
# fields.CreateRegion(region)
# fields.NodesExport("laplace_equation", "FORTRAN")
# fields.ElementsExport("laplace_equation", "FORTRAN")
# fields.Finalise()

Let the library know that you are done with computations and the resources allocated for the problem can now be released

iron.Finalise()

Visualising results

There are many different ways to visualise the results that you have generated in this example. Below we show how you can convert all the mesh data into vtk format for visualisation in paraview or other VTK friendly visualisation software.

Converting all code into vtk

import meshio
import numpy as np

# Obtain the nodeal values and their coordinates
points = np.array(NodeCoords)

# Obtain the connectivity matrix between the nodes at each element
# Note that meshio uses a 0-index node labeling system, while Opencmiss uses a 1-index node labeling system. So all node numbers need to be reduced by 1.
# elements_list[0] represents the element number
# elements_list[1:5] represents the 4 nodes associated with that element
cells_array = np.array(elements_list)[:,1:5] - 1
cells = [("tetra", cells_array)]


# Extract the values of the solution field at each of the nodes in the mesh
solution_field_nodes = np.zeros(num_nodes)
for i in range(num_nodes):
    node = i+1
    solution_field_nodes[i] = solution_field.ParameterSetGetNodeDP(
      iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,1,node,1)

# This is the code that outputs the data into a particular format
meshio.write_points_cells(
  "output_mesh.vtk",points,cells,point_data={"solution":solution_field_nodes})
    # Optionally provide extra data on points, cells, etc.
    # cell_data=cell_data,
    # field_data=field_data
# )