diff --git a/src/python/hermite_prolate_spheroid_mechanics.py b/src/python/hermite_prolate_spheroid_mechanics.py
index fec943f..a7130a9 100755
--- a/src/python/hermite_prolate_spheroid_mechanics.py
+++ b/src/python/hermite_prolate_spheroid_mechanics.py
@@ -106,7 +106,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 +125,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
@@ -279,14 +279,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..4fecef8 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:
@@ -318,11 +318,11 @@ def calcNodeAxes(i, j, k):
 
         self.componentElements = defaultdict(list)
         # Loop from endocardium to epicardium
-        for k in xrange(self.numElements[2]):
+        for k in list(range(self.numElements[2])):
             # Loop from apex to base
-            for j in xrange(self.numElements[1]):
+            for j in list(range(self.numElements[1])):
                 # Loop circumferentially
-                for i in xrange(self.numElements[0]):
+                for i in list(range(self.numElements[0])):
                     nodeAxes = calcNodeAxes(i, j, k)
                     for component in self.interpolations:
                         positions = nodeAxes[component]
@@ -360,14 +360,14 @@ def _calculateNodes(self):
         nodeNumber = itertools.count(1)
 
         # Loop from endocardium to epicardium
-        for k in xrange(self.numElements[2] * elFactor + 1):
+        for k in list(range(self.numElements[2] * elFactor + 1)):
             # Apex nodes at mu=0, have multiple apex nodes at same
             # position but different theta, get constrained to be
             # mapped together later:
             # Loop from apex to base
-            for j in xrange(self.numElements[1] * elFactor + 1):
+            for j in list(range(self.numElements[1] * elFactor + 1)):
                 # Loop circumferentially
-                for i in xrange(self.numElements[0] * elFactor):
+                for i in list(range(self.numElements[0] * elFactor)):
                     self.nodePositions.append(xyz(self.focus,
                         lambdaPositions[k], muPositions[j], thetaPositions[i]))
                     self.positionNode[(i, j, k)] = next(nodeNumber)
@@ -546,12 +546,12 @@ def constrainedNodes(self):
             elFactor = 1
         nodeSets = []
         # Loop radially
-        for k in xrange(self.numElements[2] * elFactor + 1):
+        for k in lsit(range(self.numElements[2] * elFactor + 1)):
             nodeSets.append([])
             # Longitudinal position at apex:
             j = 0
             # Loop circumferentially
-            for i in xrange(self.numElements[0] * elFactor):
+            for i in list(range(self.numElements[0] * elFactor)):
                 nodeNumber = self.nodeAtPosition((i, j, k))
                 nodeSets[-1].append(nodeNumber)
         return nodeSets