7. Using OpenCMISS from Python¶
Section author: Adam Reeve
In this tutorial we will walk through how to solve Laplace's equation on a 2D geometry using the Python bindings to OpenCMISS.
7.1. Getting Started¶
The Python code given here follows the Python Laplace example
found in the OpenCMISS examples repository. You can run the code
interactively by entering it directly into the Python interpreter,
which can be started by running python
.
Alternatively you can enter the code into a text file with a .py
extension and run it using the python command, eg:
python LaplaceExample.py
In order to use OpenCMISS we have to first import the CMISS
module from
the opencmiss
package:
from opencmiss import CMISS
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 CMISS
namespace.
7.2. Create a Coordinate System¶
First we construct a coordinate system that will be used to describe the geometry in our problem. The 2D geometry will exist in a 3D space, so we need a 3D coordinate system.
When creating an object in OpenCMISS there are at least three steps. First we initialise the object:
coordinateSystem = CMISS.CoordinateSystem()
This creates a thin wrapper that just points to the actual coordinate
system used internally by OpenCMISS, and initially the pointer is null.
Trying to use the coordinate system now would raise an exception.
To actually construct a coordinate system we call the CreateStart
method and pass it a user number.
The user number must be unique between objects of the same type and
can be used to identify the coordinate system later.
Most OpenCMISS classes require a user number when creating them,
and many also require additional parameters to the CreateStart
method:
coordinateSystemUserNumber = 1
coordinateSystem.CreateStart(coordinateSystemUserNumber)
We can now optionally set any properties on the object. We will set the dimension so that the coordinate system is 3D:
coordinateSystem.dimension = 3
And finally, we finish creating the coordinate system:
coordinateSystem.CreateFinish()
7.3. Create 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.
We can also give the region a label:
regionUserNumber = 1
region = CMISS.Region()
region.CreateStart(regionUserNumber, CMISS.WorldRegion)
region.label = "LaplaceRegion"
region.coordinateSystem = coordinateSystem
region.CreateFinish()
7.4. Create a Basis¶
The finite element description of our fields requires a basis function to interpolate field values over elements, so we create a 2D basis with linear Lagrange interpolation in both xi directions:
basisUserNumber = 1
basis = CMISS.Basis()
basis.CreateStart(basisUserNumber)
basis.type = CMISS.BasisTypes.LAGRANGE_HERMITE_TP
basis.numberOfXi = 2
basis.interpolationXi = [
CMISS.BasisInterpolationSpecifications.LINEAR_LAGRANGE,
CMISS.BasisInterpolationSpecifications.LINEAR_LAGRANGE,
]
basis.quadratureNumberOfGaussXi = [2, 2]
basis.CreateFinish()
When we set the basis type we select a value from the BasisTypes
enum.
To get a list of possible basis types and interpolation types you can
use the help
function in the python interpreter:
help(CMISS.BasisTypes)
help(CMISS.BasisInterpolationSpecifications)
Similarly, you can see what methods and properties are available for the various CMISS classes and get help information for these:
help(CMISS.Basis)
help(CMISS.Basis.CreateStart)
help(CMISS.Basis.interpolationXi)
7.5. Create a Decomposed Mesh¶
In order to define a simple 2D geometry for our problem we can use one of OpenCMISS's inbuilt generated meshes. We will create a 2D, rectangular mesh with 10 elements in both the x and y directions and tell it to use the basis we created previously:
generatedMeshUserNumber = 1
numberGlobalXElements = 10
numberGlobalYElements = 10
width = 1.0
length = 1.0
generatedMesh = CMISS.GeneratedMesh()
generatedMesh.CreateStart(generatedMeshUserNumber, region)
generatedMesh.type = CMISS.GeneratedMeshTypes.REGULAR
generatedMesh.basis = [basis]
generatedMesh.extent = [width, length, 0.0]
generatedMesh.numberOfElements = [
numberGlobalXElements,
numberGlobalYElements]
When setting the basis
property, we assign a list of
bases as we might want to construct a mesh with multiple
components using different interpolation schemes.
The generated mesh is not itself a mesh, but is used to create
a mesh. We construct the mesh
when we call the CreateFinish
method of the generated mesh
and pass in the mesh to generate:
meshUserNumber = 1
mesh = CMISS.Mesh()
generatedMesh.CreateFinish(meshUserNumber, mesh)
Here we have initialised a mesh but not called CreateStart
or CreateFinish
, instead the mesh creation is done
when finishing the creation of the generated mesh.
Because OpenCMISS can solve problems on multiple computational nodes, it must work with a decomposed mesh. We now decompose our mesh by getting the number of computational nodes and creating a decomposition with that number of domains:
decompositionUserNumber = 1
decomposition = CMISS.Decomposition()
decomposition.CreateStart(decompositionUserNumber, mesh)
decomposition.type = CMISS.DecompositionTypes.CALCULATED
decomposition.numberOfDomains = CMISS.ComputationalNumberOfNodesGet()
decomposition.CreateFinish()
Note that even when we have just one computational node, OpenCMISS still needs to work with a decomposed mesh, which will have one domain.
7.6. Defining Geometry¶
Now that we have a decomposed mesh, we can begin defining the fields we need on it. First we will create a geometric field to define our problem geometry:
geometricFieldUserNumber = 1
geometricField = CMISS.Field()
geometricField.CreateStart(geometricFieldUserNumber, region)
geometricField.meshDecomposition = decomposition
geometricField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.U, 1, 1)
geometricField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.U, 2, 1)
geometricField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.U, 3, 1)
geometricField.CreateFinish()
The call to the ComponentMeshComponentSet
method is not actually required
here as all field components will default to use the first mesh component, but
if we have defined a mesh that has multiple components (that use different interpolation
schemes) then different field components can use different mesh components.
For example, in a finite elasticity problem we could define our geometry using
quadratic Lagrange interpolation, and the hydrostatic pressure using linear Lagrange
interpolation.
We have created a field but all the field component values are currently set to zero. We can define the geometry using the generated mesh we created earlier:
generatedMesh.GeometricParametersCalculate(geometricField)
7.7. Setting up Equations¶
Now we have a geometric field we can construct an
equations set.
This defines the set of equations that we wish to solve in our
problem on this region.
The specific equation set we are solving is defined by
the fourth, fifth and sixth parameters to the CreateStart
method. These are the equations set class, type and subtype
respectively. In this example we are solving the standard Laplace
equation which is a member of the classical field equations set
class and the Laplace equation type.
When we create an equations set we also have to create an
equations set field, however, this is only used to identify
multiple equations sets of the same type on a region
so we will not use it:
equationsSetUserNumber = 1
equationsSetFieldUserNumber = 2
equationsSetField = CMISS.Field()
equationsSet = CMISS.EquationsSet()
equationsSet.CreateStart(equationsSetUserNumber, region, geometricField,
CMISS.EquationsSetClasses.CLASSICAL_FIELD,
CMISS.EquationsSetTypes.LAPLACE_EQUATION,
CMISS.EquationsSetSubtypes.STANDARD_LAPLACE,
equationsSetFieldUserNumber, equationsSetField)
equationsSet.CreateFinish()
Now we use our equations set to create a dependent field. This stores the solution to our equations:
dependentFieldUserNumber = 3
dependentField = CMISS.Field()
equationsSet.DependentCreateStart(dependentFieldUserNumber, dependentField)
equationsSet.DependentCreateFinish()
We haven't used the Field.CreateStart
method to construct
the dependent field but have had it automatically constructed by
the equations set.
We can initialise our solution with a value we think will
be close to the final solution. A field in OpenCMISS can contain multiple
field variables,
and each field variable can have multiple
components.
For the standard Laplace equation, the dependent field only has a U
variable
which has one component. Field variables can also have different field
parameter sets,
for example we can store values at a previous time step in dynamic problems.
In this example we are only interested in the VALUES
parameter set:
componentNumber = 1
initialValue = 0.5
dependentField.ComponentValuesInitialiseDP(
CMISS.FieldVariableTypes.U,
CMISS.FieldParameterSetTypes.VALUES,
componentNumber, initialValue)
Once the equations set is defined, we create the equations that use our fields to construct equations matrices and vectors. We will use sparse matrices to store the equations and enable matrix output when assembling the equations:
equations = CMISS.Equations()
equationsSet.EquationsCreateStart(equations)
equations.sparsityType = CMISS.EquationsSparsityTypes.SPARSE
equations.outputType = CMISS.EquationsOutputTypes.MATRIX
equationsSet.EquationsCreateFinish()
7.8. Defining the Problem¶
Now that we have defined all the equations we will need we can create our problem to solve. We create a standard Laplace problem, which is a member of the classical field problem class and Laplace equation problem type:
problemUserNumber = 1
problem = CMISS.Problem()
problem.CreateStart(problemUserNumber)
problem.SpecificationSet(CMISS.ProblemClasses.CLASSICAL_FIELD,
CMISS.ProblemTypes.LAPLACE_EQUATION,
CMISS.ProblemSubTypes.STANDARD_LAPLACE)
problem.CreateFinish()
The problem type defines a control loop
structure that is used when solving the problem.
We may have multiple control loops with nested sub loops,
and control loops can have different types,
for example load incremented loops or time loops for dynamic problems.
In this example a simple, single iteration loop is created without any sub loops.
If we wanted to access the control loop and modify
it we would use the problem.ControlLoopGet
method before
finishing the creation of the control loops, but we will just
leave it with the default configuration:
problem.ControlLoopCreateStart()
problem.ControlLoopCreateFinish()
7.9. Configuring Solvers¶
After defining the problem structure we can create the
solvers that
will be run to actually solve our problem.
The problem type defines the solvers to be set up
so we call problem.SolversCreatStart
to create the solvers and
then we can access the solvers to modify their properties:
solver = CMISS.Solver()
problem.SolversCreateStart()
problem.SolverGet([CMISS.ControlLoopIdentifiers.NODE], 1, solver)
solver.outputType = CMISS.SolverOutputTypes.SOLVER
solver.linearType = CMISS.LinearSolverTypes.ITERATIVE
solver.linearIterativeAbsoluteTolerance = 1.0e-10
solver.linearIterativeRelativeTolerance = 1.0e-10
problem.SolversCreateFinish()
Note that we initialised a solver but didn't create it directly
by calling its CreateStart
method,
it was created with the call to SolversCreateStart
and then we obtain
it with the call to SolverGet
. If we look at the help for the
SolverGet
method we see it takes three parameters:
- controlLoopIdentifiers
- A list of integers used to identify the control loop to get a solver for.
This always starts with the root control loop, given by
CMISS.ControlLoopIdentifiers.NODE
. In this example we only have the one control loop and no sub loops. - solverIndex
- The index of the solver to get, as a control loop may have multiple solvers. In this case there is only one solver in our root control loop.
- solver
- An initialised solver object that hasn't been created yet, and on return it will be the solver that we asked for.
Once we've obtained the solver we then set various properties before finishing the creation of all the problem solvers.
After defining our solver we can create the equations for the solver to solve by adding our equations sets to the solver equations. In this example we have just one equations set to add but for coupled problems we may have multiple equations sets in the solver equations. We also tell OpenCMISS to use sparse matrices to store our solver equations:
solverEquations = CMISS.SolverEquations()
problem.SolverEquationsCreateStart()
solver.SolverEquationsGet(solverEquations)
solverEquations.sparsityType = CMISS.SolverEquationsSparsityTypes.SPARSE
equationsSetIndex = solverEquations.EquationsSetAdd(equationsSet)
problem.SolverEquationsCreateFinish()
7.10. Setting Boundary Conditions¶
The final step in configuring the problem is to define the boundary conditions to be satisfied. We will set the dependent field value at the first node to be zero, and at the last node to be 1.0. These nodes will correspond to opposite corners in our geometry. Because OpenCMISS can solve our problem on multiple computational nodes where each computational node does not necessarily know about all nodes in our mesh, we must first check that the node we are setting the boundary condition at is in our computational node domain:
nodes = CMISS.Nodes()
region.NodesGet(nodes)
firstNodeNumber = 1
lastNodeNumber = nodes.numberOfNodes
firstNodeDomain = decomposition.NodeDomainGet(firstNodeNumber, 1)
lastNodeDomain = decomposition.NodeDomainGet(lastNodeNumber, 1)
computationalNodeNumber = CMISS.ComputationalNodeNumberGet()
boundaryConditions = CMISS.BoundaryConditions()
solverEquations.BoundaryConditionsCreateStart(boundaryConditions)
if firstNodeDomain == computationalNodeNumber:
boundaryConditions.SetNode(
dependentField, CMISS.FieldVariableTypes.U,
1, 1, firstNodeNumber, 1,
CMISS.BoundaryConditionsTypes.FIXED, 0.0)
if lastNodeDomain == computationalNodeNumber:
boundaryConditions.SetNode(
dependentField, CMISS.FieldVariableTypes.U,
1, 1, lastNodeNumber, 1,
CMISS.BoundaryConditionsTypes.FIXED, 1.0)
solverEquations.BoundaryConditionsCreateFinish()
When setting a boundary condition at a node we can use either the AddNode
method or the SetNode
method. Using AddNode
will add the value
we provide to the current field value and set this as the boundary condition value,
but here we want to directly specify the value so we use the SetNode
method.
The arguments to the SetNode
method are the field, field variable type,
node version number, node user number, node derivative number, field component number,
boundary condition type and boundary condition value.
The version and derivative numbers are one as we aren't using versions and we are
setting field values rather than derivative values.
We can also only set derivative boundary conditions when using a Hermite basis type.
There are a wide number of boundary condition types that can be set but
many are only available for certain equation set types and in this example we
simply want to fix the field value.
When solverEquations.BoundaryConditionsCreateFinish()
is called
OpenCMISS will construct the solver matrices and vectors.
7.11. Solving¶
After our problem solver equations have been fully defined we are now ready
to solve our problem. When we call the Solve
method of the problem it will
loop over the control loops and control loop solvers to solve our problem:
problem.Solve()
7.12. Exporting the Solution¶
Once the problem has been solved, the dependent field contains the solution
to our problem. We can then export the dependent and geometric fields to a
FieldML file so that we can visualise the solution using cmgui.
We will export the geometric and dependent field values to
a LaplaceExample.xml
file.
Separate plain text data files will also be created:
baseName = "laplace"
dataFormat = "PLAIN_TEXT"
fml = CMISS.FieldMLIO()
fml.OutputCreate(mesh, "", baseName, dataFormat)
fml.OutputAddFieldNoType(
baseName + ".geometric", dataFormat, geometricField,
CMISS.FieldVariableTypes.U, CMISS.FieldParameterSetTypes.VALUES)
fml.OutputAddFieldNoType(
baseName + ".phi", dataFormat, dependentField,
CMISS.FieldVariableTypes.U, CMISS.FieldParameterSetTypes.VALUES)
fml.OutputWrite("LaplaceExample.xml")
fml.Finalise()