diff --git a/src/python/hermite_prolate_spheroid_mechanics.py b/src/python/hermite_prolate_spheroid_mechanics.py
index fec943f..47b36fa 100755
--- a/src/python/hermite_prolate_spheroid_mechanics.py
+++ b/src/python/hermite_prolate_spheroid_mechanics.py
@@ -59,21 +59,26 @@
 equationsSetUserNumber = 1
 problemUserNumber = 1
 
+worldRegion = iron.Region()
+iron.Context.WorldRegionGet(worldRegion)
+
 # Get the number of computational nodes and this computational node number
 # for when running in parallel with MPI
-numberOfComputationalNodes = iron.ComputationalNumberOfNodesGet()
-computationalNodeNumber = iron.ComputationalNodeNumberGet()
+computationEnvironment = iron.ComputationEnvironment()
+iron.Context.ComputationEnvironmentGet(computationEnvironment)
+numberOfComputationalNodes = computationEnvironment.NumberOfWorldNodesGet()
+computationalNodeNumber = computationEnvironment.WorldNodeNumberGet()
 
 # Create a 3D rectangular cartesian coordinate system
 coordinateSystem = iron.CoordinateSystem()
-coordinateSystem.CreateStart(coordinateSystemUserNumber)
+coordinateSystem.CreateStart(coordinateSystemUserNumber,iron.Context)
 coordinateSystem.DimensionSet(3)
 coordinateSystem.CreateFinish()
 
 # Create a region within the world region and
 # assign the coordinate system to the region
 region = iron.Region()
-region.CreateStart(regionUserNumber, iron.WorldRegion)
+region.CreateStart(regionUserNumber, worldRegion)
 region.LabelSet("ProlateSpheroid")
 region.CoordinateSystemSet(coordinateSystem)
 region.CreateFinish()
@@ -106,7 +111,7 @@
 geometricField.CreateFinish()
 
 # Use the prolate spheroid geometry to set the geometric field parameters
-geometry.setGeometry(geometricField)
+geometry.setGeometry(computationEnvironment,geometricField)
 
 # Create a fibre field and attach it to the geometric field
 # This has three components describing fibre orientations as
@@ -125,7 +130,7 @@
 fibreField.CreateFinish()
 
 # Use the prolate spheroid geometry to set up the fibre field values
-geometry.setFibres(fibreField)
+geometry.setFibres(computationEnvironment,fibreField)
 
 # Create the equations set and equations set field
 # This defines the type of equations to solve
@@ -231,7 +236,7 @@
 problemSpecification = [iron.ProblemClasses.ELASTICITY,
         iron.ProblemTypes.FINITE_ELASTICITY,
         iron.ProblemSubtypes.NONE]
-problem.CreateStart(problemUserNumber, problemSpecification)
+problem.CreateStart(problemUserNumber,iron.Context,problemSpecification)
 problem.CreateFinish()
 
 # Create the problem control loops
@@ -279,14 +284,14 @@
 boundaryConditions = iron.BoundaryConditions()
 solverEquations.BoundaryConditionsCreateStart(boundaryConditions)
 
-def getDomainNodes(geometry, decomposition, component):
+def getDomainNodes(computationEnvironment, geometry, decomposition, component):
     component_name = interpolations[component - 1]
-    computationalNodeNumber = iron.ComputationalNodeNumberGet()
+    computationalNodeNumber = computationEnvironment.WorldNodeNumberGet()
     nodes = geometry.componentNodes(component_name)
     meshComponent = geometry.meshComponent(component_name)
     return set(node for node in nodes
         if decomposition.NodeDomainGet(node, meshComponent) == computationalNodeNumber)
-geometricDomainNodes = getDomainNodes(geometry, decomposition, geometricMeshComponent)
+geometricDomainNodes = getDomainNodes(computationEnvironment, geometry, decomposition, geometricMeshComponent)
 
 # Fix epicardium nodes at the base:
 baseNodes = set(geometry.nodeGroup('base'))
diff --git a/src/python/prolate_spheroid_geometry.py b/src/python/prolate_spheroid_geometry.py
index af62c28..229a40e 100644
--- a/src/python/prolate_spheroid_geometry.py
+++ b/src/python/prolate_spheroid_geometry.py
@@ -136,11 +136,11 @@ def generateMesh(self, region):
         mesh.CreateFinish()
         return mesh
 
-    def setGeometry(self, geometricField):
+    def setGeometry(self, computationEnvironment, geometricField):
         decomposition = iron.Decomposition()
         geometricField.MeshDecompositionGet(decomposition)
         geometricInterpolation = self.interpolations[0]
-        compNodeNumber = iron.ComputationalNodeNumberGet()
+        compNodeNumber = computationEnvironment.WorldNodeNumberGet()
         geometricMeshComponent = 1
         # Set the geometric field parameters from the prolate spheroid geometry
         for nodeNum, values in enumerate(self.nodeValues(), 1):
@@ -169,12 +169,12 @@ def setGeometry(self, geometricField):
                             iron.FieldVariableTypes.U,
                             iron.FieldParameterSetTypes.VALUES)
 
-    def setFibres(self, fibreField):
+    def setFibres(self, computationEnvironment,fibreField):
         endocardiumFibreAngle = self.endocardiumFibreAngle
         epicardiumFibreAngle = self.epicardiumFibreAngle
         sheetAngle = self.sheetAngle
         geometricMeshComponent = 1
-        compNodeNumber = iron.ComputationalNodeNumberGet()
+        compNodeNumber = computationEnvironment.WorldNodeNumberGet()
 
         hasQuadratic = self.interpolations[0] == 'quadratic'
         if hasQuadratic:
@@ -252,7 +252,7 @@ def makeBasis(userNumber, interpolationName, collapsed):
                     }
             interpolation = cmissInterpolations[interpolationName]
             basis = iron.Basis()
-            basis.CreateStart(basisUserNumber)
+            basis.CreateStart(basisUserNumber,iron.Context)
             basis.TypeSet(iron.BasisTypes.LAGRANGE_HERMITE_TP)
             basis.NumberOfXiSet(3)
             basis.InterpolationXiSet([interpolation] * 3)