diff --git a/README.rst b/README.rst index fe1b24e..b7d999d 100644 --- a/README.rst +++ b/README.rst @@ -10,8 +10,10 @@ Building the example Instructions on how to configure and build with CMake:: git clone https://github.com/OpenCMISS-Examples/left_ventricle_inflation.git + cd left_ventricle_inflation mkdir build - cmake -DOpenCMISSLibs_DIR=/path/to/opencmisslib/install ../left_ventricle_inflation + cd build + cmake -DOpenCMISS_INSTALL_ROOT=/path/to/opencmiss/install ../. make # cmake --build . will also work here and is much more platform agnostic. Running the example @@ -19,7 +21,6 @@ Running the example Explain how the example is run:: - cd build ./src/fortran/left_ventricle_inflation.F90 or maybe it is a Python only example:: diff --git a/src/python/exfile.py b/src/python/exfile.py index 77f58f7..ab228cf 100644 --- a/src/python/exfile.py +++ b/src/python/exfile.py @@ -127,7 +127,7 @@ def _read_element(self, f): try: indices = map(int, element_line.split(':')[1].split()) except: - print element_line + print(element_line) raise if indices[1] == 0 and indices[2] == 0: #raise ExfileError(f, "Face or line elements not supported") @@ -292,12 +292,13 @@ def _read_node(self, f): while read < self.num_node_values: line = f.readline() try: - new_values = map(float, line.split()) + new_values = list(map(float, line.split())) except ValueError: raise ExfileError(f, "Expecting node values, got: %s" % line.strip()) - if read + len(new_values) > self.num_node_values: + if read + len(list(new_values)) > self.num_node_values: raise ExfileError(f, "Got more node values than expected.") - values[read:read + len(new_values)] = new_values + for nodeIndex in (0,len(new_values)-1): + values[read+nodeIndex] = new_values[nodeIndex] read += len(new_values) self.nodes.append(ExnodeNode(number, values)) @@ -462,7 +463,7 @@ def _read_element(self, f): element_line = f.readline() if element_line == "": raise EOFError - indices = map(int, element_line.split(':')[1].split()) + indices = list(map(int, element_line.split(':')[1].split())) if indices[1] == 0 and indices[2] == 0: # raise ExfileError(f, "Face or line elements not supported") values = [] diff --git a/src/python/left_ventricle_inflation.py b/src/python/left_ventricle_inflation.py index f137e95..12025de 100644 --- a/src/python/left_ventricle_inflation.py +++ b/src/python/left_ventricle_inflation.py @@ -26,7 +26,7 @@ #> of Oxford are Copyright (C) 2007 by the University of Auckland and #> the University of Oxford. All Rights Reserved. #> -#> Contributor(s): Sander Land +#> Contributor(s): Zhinuo Jenny Wang, Sander Land #> #> Alternatively, the contents of this file may be used under the terms of #> either the GNU General Public License Version 2 or later (the "GPL"), or @@ -77,13 +77,15 @@ cellMLOption = [False] ### Set arbitrary user numbers which are unique to each object. -(coordinateSystemUserNumber, +(contextUserNumber, + coordinateSystemUserNumber, regionUserNumber, linearBasisUserNumber, cubicBasisUserNumber, generatedMeshUserNumber, meshUserNumber, decompositionUserNumber, + decomposerUserNumber, geometricFieldUserNumber, equationsSetUserNumber, equationsSetFieldUserNumber, @@ -96,7 +98,7 @@ cellMLUserNumber, cellMLModelsFieldUserNumber, cellMLParametersFieldUserNumber, - cellMLIntermediateFieldUserNumber) = range(1, 21) + cellMLIntermediateFieldUserNumber) = range(1, 23) # Set up geometric model # longitudinal elements, circumferential elements, transmural elements @@ -116,23 +118,30 @@ apex_elems.append([start, end]) # Set up region and CS -[numOfCompNodes, compNodeNum, CS, region] = BasicSetUp(regionUserNumber, coordinateSystemUserNumber) +[context, worldWorkGroup, numOfCompNodes, compNodeNum, CS, region] = BasicSetUp(contextUserNumber, regionUserNumber, coordinateSystemUserNumber) + +# Set the OpenCMISS random seed so that we can test this example by using the +# same parallel decomposition +numberOfRandomSeeds = context.RandomSeedsSizeGet() +randomSeeds = [0]*numberOfRandomSeeds +randomSeeds[0] = 100 +context.RandomSeedsSet(randomSeeds) # Set up tricubic Hermite basis functions -[linearBasis, colBasis] = BasisFunction(linearBasisUserNumber, numOfXi, option, collapsed=True) +[linearBasis, colBasis] = BasisFunction(context, linearBasisUserNumber, numOfXi, option, collapsed=True) # Set up mesh -mesh = iron.Mesh() +mesh = oc.Mesh() mesh.CreateStart(meshUserNumber, region, numOfXi) mesh.NumberOfComponentsSet(1) mesh.NumberOfElementsSet(inputElems.num_elements) -nodes = iron.Nodes() +nodes = oc.Nodes() nodes.CreateStart(region, inputNodes.num_nodes) nodes.CreateFinish() # Linear lagrange component -linearElem = iron.MeshElements() +linearElem = oc.MeshElements() linearElem.CreateStart(mesh, 1, linearBasis) for elem in inputElems.elements: for i in range(1, elems[2]+1): @@ -140,7 +149,7 @@ if (elem.number >= apex_elems[i-1][0]) & (elem.number <= apex_elems[i-1][1]): linearElem.BasisSet(elem.number, colBasis) nodes = list(OrderedDict.fromkeys(elem.nodes)) - nodes = map(int, nodes) + nodes = list(map(int, nodes)) linearElem.NodesSet(elem.number, nodes) else: linearElem.NodesSet(elem.number, elem.nodes) @@ -149,13 +158,16 @@ mesh.CreateFinish() # Set up decomposition for the mesh. -decomposition = DecompositionSetUp(decompositionUserNumber, mesh, numOfCompNodes) +decomposition = DecompositionSetUp(decompositionUserNumber, mesh) + +# Decompose +decomposer = DecomposerSetUp(decomposerUserNumber, region, worldWorkGroup, decomposition) # Set up geometric field. geometricField = GeometricFieldSetUp(geometricFieldUserNumber, region, decomposition, option) # Update the geometric field parameters manually -geometricField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES) +geometricField.ParameterSetUpdateStart(oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES) def ExtractNodeCoords(nodes, field_name): coords = [] @@ -178,7 +190,7 @@ def ExtractNodeCoords(nodes, field_name): coord = [] for component in [1,2,3]: value = all_nodes[node_num-1][component-1] - geometricField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, 1, + geometricField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, 1, node_num, component, value) coord = all_nodes[node_num-1] if coord[2] >= 5: @@ -222,7 +234,7 @@ def ExtractNodeCoords(nodes, field_name): eval_node_num.append(i+1) #print eval_node_num -geometricField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES) +geometricField.ParameterSetUpdateFinish(oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES) # Export undeformed geometry. GeometricFieldExport(region, "LVInflation_trilinear_undeformed_"+str(elems[2])+"-"+str(elems[1])+"-"+str(elems[0])) @@ -237,15 +249,15 @@ def ExtractNodeCoords(nodes, field_name): cellMLOption) # Set up equations set -equationsSetField = iron.Field() # Equations are also in a field -equationsSet = iron.EquationsSet() # Initialise an equation set. -equationsSetSpecification = [iron.EquationsSetClasses.ELASTICITY, - iron.EquationsSetTypes.FINITE_ELASTICITY, - iron.EquationsSetSubtypes.TRANSVERSE_ISOTROPIC_GUCCIONE] +equationsSetField = oc.Field() # Equations are also in a field +equationsSet = oc.EquationsSet() # Initialise an equation set. +equationsSetSpecification = [oc.EquationsSetClasses.ELASTICITY, + oc.EquationsSetTypes.FINITE_ELASTICITY, + oc.EquationsSetSubtypes.TRANSVERSE_ISOTROPIC_GUCCIONE] equationsSet.CreateStart(equationsSetUserNumber, region, fibreField, equationsSetSpecification, equationsSetFieldUserNumber, equationsSetField) equationsSet.CreateFinish() -print "----> Set up equations set <---\n" +print("----> Set up equations set <---\n") # Set up material field in equations set. equationsSet.MaterialsCreateStart(materialFieldUserNumber, materialField) @@ -270,9 +282,9 @@ def ExtractNodeCoords(nodes, field_name): increm = pressure_increments[i] p = p + increm tol = tolerances[i] - print 'Applying pressure increment of: ', increm, ' using ', iters, ' iterations' - print 'Current pressure is: ', p - [problem, solverEquations] = ProblemSolverSetup(equationsSet, problemUserNumber, iters, tol, cellMLOption) + print('Applying pressure increment of: ', increm, ' using ', iters, ' iterations') + print('Current pressure is: ', p) + [problem, solverEquations] = ProblemSolverSetup(context, equationsSet, problemUserNumber, iters, tol, cellMLOption) BCEndoPressure(solverEquations, dependentField, endocardial_nodes, increm, basal_nodes, option) # Solve Problem diff --git a/src/python/lib.py b/src/python/lib.py index e9d4f90..bd6f49a 100644 --- a/src/python/lib.py +++ b/src/python/lib.py @@ -12,90 +12,100 @@ # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# -from opencmiss.iron import iron +from opencmiss.opencmiss import OpenCMISS_Python as oc import numpy import math import os # =================================================================================# -def BasicSetUp(regionUserNumber, coordinateSystemUserNumber): - # This function sets up the world region, 3D CS, parallel computing nodes, and - # diagnostics. +def BasicSetUp(contextUserNumber, regionUserNumber, coordinateSystemUserNumber): + # This function sets up the context, world region, 3D CS, parallel computing + # nodes, and diagnostics. + + context = oc.Context() + context.Create(contextUserNumber) # Set up diagnostics/debug - #iron.DiagnosticsSetOn(iron.DiagnosticTypes.IN,[1,2,3,4,5], + #oc.DiagnosticsSetOn(oc.DiagnosticTypes.IN,[1,2,3,4,5], #"Diagnostics",["DOMAIN_MAPPINGS_LOCAL_FROM_GLOBAL_CALCULATE"]) # Get computational node information for parallel computing - numberOfComputationalNodes = iron.ComputationalNumberOfNodesGet() - computationalNodeNumber = iron.ComputationalNodeNumberGet() + computationEnvironment = oc.ComputationEnvironment() + context.ComputationEnvironmentGet(computationEnvironment) + worldWorkGroup = oc.WorkGroup() + computationEnvironment.WorldWorkGroupGet(worldWorkGroup) + numberOfComputationalNodes = worldWorkGroup.NumberOfGroupNodesGet() + computationalNodeNumber = worldWorkGroup.GroupNodeNumberGet() # Set up 3D RC coordinate system - coordinateSystem = iron.CoordinateSystem() - coordinateSystem.CreateStart(coordinateSystemUserNumber) + coordinateSystem = oc.CoordinateSystem() + coordinateSystem.CreateStart(coordinateSystemUserNumber,context) coordinateSystem.dimension = 3 coordinateSystem.CreateFinish() + worldRegion = oc.Region() + context.WorldRegionGet(worldRegion) + # Create world region - region = iron.Region() - region.CreateStart(regionUserNumber, iron.WorldRegion) + region = oc.Region() + region.CreateStart(regionUserNumber, worldRegion) region.label = "Region" region.coordinateSystem = coordinateSystem region.CreateFinish() # Output for diagnostics - print "----> Set up coordinate system and world region <----\n" + print("----> Set up coordinate system and world region <----\n") - return numberOfComputationalNodes, computationalNodeNumber, coordinateSystem, region + return context, worldWorkGroup, numberOfComputationalNodes, computationalNodeNumber, coordinateSystem, region # =================================================================================# #=================================================================================# -def BasisFunction(basisUserNumber, numOfXi, option, collapsed): +def BasisFunction(context, basisUserNumber, numOfXi, option, collapsed): # This function sets up the basis function depending on the option given. if option[0] == 1: # Trilinear basis function for interpolation of geometry. - basis = iron.Basis() - basis.CreateStart(basisUserNumber) + basis = oc.Basis() + basis.CreateStart(basisUserNumber,context) basis.numberOfXi = numOfXi - basis.type = iron.BasisTypes.LAGRANGE_HERMITE_TP - basis.interpolationXi = [iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE] * numOfXi + basis.type = oc.BasisTypes.LAGRANGE_HERMITE_TP + basis.interpolationXi = [oc.BasisInterpolationSpecifications.LINEAR_LAGRANGE] * numOfXi basis.QuadratureLocalFaceGaussEvaluateSet(True) basis.quadratureNumberOfGaussXi = [2,2,2] basis.CreateFinish() # Output for diagnostics - print "----> Set up trilinear basis functions for geometry, use element based interpolation for pressure <----\n" + print("----> Set up trilinear basis functions for geometry, use element based interpolation for pressure <----\n") if collapsed: - basisCol = iron.Basis() - basisCol.CreateStart(basisUserNumber+1) + basisCol = oc.Basis() + basisCol.CreateStart(basisUserNumber+1,context) basisCol.numberOfXi = numOfXi - basisCol.type = iron.BasisTypes.LAGRANGE_HERMITE_TP - basisCol.interpolationXi = [iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE] * numOfXi + basisCol.type = oc.BasisTypes.LAGRANGE_HERMITE_TP + basisCol.interpolationXi = [oc.BasisInterpolationSpecifications.LINEAR_LAGRANGE] * numOfXi basisCol.QuadratureLocalFaceGaussEvaluateSet(True) basisCol.quadratureNumberOfGaussXi = [2,2,2] - basisCol.CollapsedXiSet([iron.BasisXiCollapse.XI_COLLAPSED, iron.BasisXiCollapse.COLLAPSED_AT_XI0, iron.BasisXiCollapse.NOT_COLLAPSED]) - print "---> Set up collapsed basis functions for apical elements" + basisCol.CollapsedXiSet([oc.BasisXiCollapse.XI_COLLAPSED, oc.BasisXiCollapse.COLLAPSED_AT_XI0, oc.BasisXiCollapse.NOT_COLLAPSED]) + print("---> Set up collapsed basis functions for apical elements") basisCol.CreateFinish() return basis, basisCol return basis elif option[0] == 2: - quadBasis = iron.Basis() - quadBasis.CreateStart(basisUserNumber[0]) - quadBasis.InterpolationXiSet([iron.BasisInterpolationSpecifications.QUADRATIC_LAGRANGE]*numOfXi) + quadBasis = oc.Basis() + quadBasis.CreateStart(basisUserNumber[0],context) + quadBasis.InterpolationXiSet([oc.BasisInterpolationSpecifications.QUADRATIC_LAGRANGE]*numOfXi) quadBasis.QuadratureNumberOfGaussXiSet([4]*numOfXi) quadBasis.QuadratureLocalFaceGaussEvaluateSet(True) quadBasis.CreateFinish() # Tricubic Hermite basis function for interpolation of geometry. - cubicBasis = iron.Basis() # For geometry. - cubicBasis.CreateStart(basisUserNumber[1]) - cubicBasis.InterpolationXiSet([iron.BasisInterpolationSpecifications.CUBIC_HERMITE] * numOfXi) + cubicBasis = oc.Basis() # For geometry. + cubicBasis.CreateStart(basisUserNumber[1],context) + cubicBasis.InterpolationXiSet([oc.BasisInterpolationSpecifications.CUBIC_HERMITE] * numOfXi) cubicBasis.QuadratureNumberOfGaussXiSet([4] * numOfXi) cubicBasis.QuadratureLocalFaceGaussEvaluateSet(True) cubicBasis.CreateFinish() # Output for diagnostics - print "----> Set up tricubic hermite basis function for geometry and trilinear for hydrostatic pressure <----\n" + print("----> Set up tricubic hermite basis function for geometry and trilinear for hydrostatic pressure <----\n") return quadBasis, cubicBasis @@ -104,13 +114,13 @@ def BasisFunction(basisUserNumber, numOfXi, option, collapsed): #=================================================================================# def GeneratedMesh(generatedMeshUserNumber, meshUserNumber, region, bases, dimensions, elements): # This function sets up a generated mesh using user specified dimensions. - generatedMesh = iron.GeneratedMesh() + generatedMesh = oc.GeneratedMesh() generatedMesh.CreateStart(generatedMeshUserNumber, region) - generatedMesh.TypeSet(iron.GeneratedMeshTypes.REGULAR) + generatedMesh.TypeSet(oc.GeneratedMeshTypes.REGULAR) generatedMesh.BasisSet(bases) generatedMesh.ExtentSet(dimensions) generatedMesh.NumberOfElementsSet(elements) - mesh = iron.Mesh() + mesh = oc.Mesh() generatedMesh.CreateFinish(meshUserNumber, mesh) return generatedMesh, mesh @@ -119,40 +129,54 @@ def GeneratedMesh(generatedMeshUserNumber, meshUserNumber, region, bases, dimens #=================================================================================# #=================================================================================# -def DecompositionSetUp(decompositionUserNumber, mesh, numberOfComputationalNodes): +def DecompositionSetUp(decompositionUserNumber, mesh): # This function sets up the decomposition of the mesh. - decomposition = iron.Decomposition() + decomposition = oc.Decomposition() decomposition.CreateStart(decompositionUserNumber, mesh) - decomposition.type = iron.DecompositionTypes.CALCULATED - decomposition.NumberOfDomainsSet(numberOfComputationalNodes) + decomposition.type = oc.DecompositionTypes.CALCULATED decomposition.CalculateFacesSet(True) decomposition.CreateFinish() # Output for diagnostics - print "----> Set up decomposition <----\n" + print("----> Set up decomposition <----\n") return decomposition +#=================================================================================# + +#=================================================================================# +def DecomposerSetUp(decomposerUserNumber, region, workGroup, decomposition): + # This function sets up the decomposer. + decomposer = oc.Decomposer() + decomposer.CreateStart(decomposerUserNumber, region, workGroup) + decompositionIndex = decomposer.DecompositionAdd(decomposition) + decomposer.CreateFinish() + + # Output for diagnostics + print("----> Set up decomposer <----\n") + return decomposer + + #=================================================================================# #=================================================================================# def GeometricFieldSetUp(geometricFieldUserNumber, region, decomposition, option): # Set up geometry field - geometricField = iron.Field() # Initialise + geometricField = oc.Field() # Initialise geometricField.CreateStart(geometricFieldUserNumber, region) - geometricField.MeshDecompositionSet(decomposition) - geometricField.VariableLabelSet(iron.FieldVariableTypes.U, "Geometry") + geometricField.DecompositionSet(decomposition) + geometricField.VariableLabelSet(oc.FieldVariableTypes.U, "Geometry") if option[0] == 2: # Tricubic Hermite if option[1] == 1: - geometricField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + geometricField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) # Output for diagnostics - print "----> Set up tricubic Hermite geometric field with unit scaling <----\n" + print("----> Set up tricubic Hermite geometric field with unit scaling <----\n") elif option[1] == 2: - geometricField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + geometricField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) # Output for diagnostics - print "----> Set up tricubic Hermite geometric field with arithmetic mean scaling <----\n" + print("----> Set up tricubic Hermite geometric field with arithmetic mean scaling <----\n") geometricField.CreateFinish() @@ -166,28 +190,28 @@ def GeometricFieldInitialise(xNodes, yNodes, zNodes, geometricField, numNodes, o # This function initialises the geometric field with user specified coordinates. # Initialise nodal values. for node, value in enumerate(xNodes, 1): - geometricField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 1, value) + geometricField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 1, value) for node, value in enumerate(yNodes, 1): - geometricField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 2, value) + geometricField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 2, value) for node, value in enumerate(zNodes, 1): - geometricField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 3, value) + geometricField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 3, value) # Initialise first derivatives. if option[0] == 2: # Tricubic Hermite basis. for node in range(numNodes): - geometricField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1, node + 1, 1, max(xNodes)) - geometricField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2, node + 1, 2, max(yNodes)) - geometricField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3, node + 1, 3, max(zNodes)) + geometricField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1, node + 1, 1, max(xNodes)) + geometricField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2, node + 1, 2, max(yNodes)) + geometricField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S3, node + 1, 3, max(zNodes)) # Output - print "----> Initialised geometric nodal values <----\n" + print("----> Initialised geometric nodal values <----\n") return geometricField @@ -201,14 +225,14 @@ def GeometricFieldExport(region, filename): if not os.path.exists("./results"): os.makedirs("./results") - exportField = iron.Fields() + exportField = oc.Fields() exportField.CreateRegion(region) exportField.NodesExport("./results/" + filename, "FORTRAN") exportField.ElementsExport("./results/" + filename, "FORTRAN") exportField.Finalise() # Output - print "----> Export undeformed geometry <----\n" + print("----> Export undeformed geometry <----\n") #=================================================================================# @@ -219,13 +243,13 @@ def ExtractNodesElements(filename): try: fid_node = open(filename+'.exnode', 'r') except IOError: - print 'ERROR: Unable to open '+filename+'.exnode' + print('ERROR: Unable to open '+filename+'.exnode') return try: fid_elem = open(filename+'.exelem', 'r') except IOError: - print 'ERROR: Unable to open '+filename+'.exelem' + print('ERROR: Unable to open '+filename+'.exelem') return for i in range(1,86): @@ -279,45 +303,45 @@ def ExtractNodesElements(filename): def FibreFieldSetUp(fibreFieldUserNumber, region, decomposition, geometricField, option, microstructure, inputNodes): # This function sets up the fibre field and initialises the values. # Sets up the fibre field. - fibreField = iron.Field() + fibreField = oc.Field() fibreField.CreateStart(fibreFieldUserNumber, region) - fibreField.TypeSet(iron.FieldTypes.FIBRE) - fibreField.MeshDecompositionSet(decomposition) + fibreField.TypeSet(oc.FieldTypes.FIBRE) + fibreField.DecompositionSet(decomposition) fibreField.GeometricFieldSet(geometricField) - fibreField.VariableLabelSet(iron.FieldVariableTypes.U, "Fibre") + fibreField.VariableLabelSet(oc.FieldVariableTypes.U, "Fibre") if option[0] == 1: fibreField.NumberOfVariablesSet(1) - fibreField.NumberOfComponentsSet(iron.FieldVariableTypes.U, 3) + fibreField.NumberOfComponentsSet(oc.FieldVariableTypes.U, 3) if microstructure == 1: for component in [1, 2, 3]: - fibreField.ComponentInterpolationSet(iron.FieldVariableTypes.U, component, - iron.FieldInterpolationTypes.CONSTANT) + fibreField.ComponentInterpolationSet(oc.FieldVariableTypes.U, component, + oc.FieldInterpolationTypes.CONSTANT) elif microstructure == 2: for component in [1, 2, 3]: - fibreField.ComponentMeshComponentSet(iron.FieldVariableTypes.U, component, 1) + fibreField.ComponentMeshComponentSet(oc.FieldVariableTypes.U, component, 1) elif option[0] == 2: # Tricubic Hermite interpolation if option[1] == 1: - fibreField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + fibreField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) # Output - print "----> Set up tricubic hermite fibre field with unit scaling <----\n" + print("----> Set up tricubic hermite fibre field with unit scaling <----\n") elif option[1] == 2: - fibreField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + fibreField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) # Output - print "----> Set up tricubic hermite fibre field with arithmetic mean scaling <----\n" + print("----> Set up tricubic hermite fibre field with arithmetic mean scaling <----\n") if microstructure == 1: # Homogeneous fibre field. for component in [1, 2, 3]: - fibreField.ComponentInterpolationSet(iron.FieldVariableTypes.U, component, - iron.FieldInterpolationTypes.CONSTANT) + fibreField.ComponentInterpolationSet(oc.FieldVariableTypes.U, component, + oc.FieldInterpolationTypes.CONSTANT) elif microstructure == 2: # Heterogeneous fibre field using linear interpolation. for component in [1, 2, 3]: - fibreField.ComponentMeshComponentSet(iron.FieldVariableTypes.U, component, 1) + fibreField.ComponentMeshComponentSet(oc.FieldVariableTypes.U, component, 1) fibreField.CreateFinish() @@ -330,10 +354,10 @@ def FibreFieldSetUp(fibreFieldUserNumber, region, decomposition, geometricField, angle = inputNodes.node_values("fibers", component_name, n) angle = float(angle[0]) angle = angle*math.pi/180 - fibreField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, - 1, iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, n, + fibreField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, + 1, oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, n, component, angle) - print "----> Initialised heterogeneous fibre angles <----\n" + print("----> Initialised heterogeneous fibre angles <----\n") return fibreField @@ -344,24 +368,24 @@ def MaterialFieldSetUpAuto(materialFieldUserNumber, equationsSet, params, cellML # This function is used for setting up material field when using CellML # description of constitutive model. # Sets up material field, and apply field to mesh component. - materialField = iron.Field() + materialField = oc.Field() equationsSet.MaterialsCreateStart(materialFieldUserNumber, materialField) - materialField.VariableLabelSet(iron.FieldVariableTypes.U, "Material") + materialField.VariableLabelSet(oc.FieldVariableTypes.U, "Material") if cellMLOption[0]: - print "----> CellML Material Field using gauss point interpolation <----\n" + print("----> CellML Material Field using gauss point interpolation <----\n") for component, param in enumerate(params, 1): - materialField.ComponentInterpolationSet(iron.FieldVariableTypes.U, component, - iron.FieldInterpolationTypes.GAUSS_POINT_BASED) + materialField.ComponentInterpolationSet(oc.FieldVariableTypes.U, component, + oc.FieldInterpolationTypes.GAUSS_POINT_BASED) materialField.CreateFinish() ######################################################################### # Initialise parameter values. for component, param in enumerate(params, 1): - materialField.ComponentValuesInitialiseDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, + materialField.ComponentValuesInitialiseDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, component, param) # Output - print "----> Initialised " + str(len(params)) + " material parameters <----\n" + print("----> Initialised " + str(len(params)) + " material parameters <----\n") return materialField, equationsSet @@ -370,44 +394,44 @@ def MaterialFieldSetUpAuto(materialFieldUserNumber, equationsSet, params, cellML #=================================================================================# def MaterialFieldSetUp(materialFieldUserNumber, region, decomposition, geometricField, params, option, cellMLOption): # Sets up material field, and apply field to mesh component. - materialField = iron.Field() + materialField = oc.Field() materialField.CreateStart(materialFieldUserNumber, region) - materialField.TypeSet(iron.FieldTypes.MATERIAL) - materialField.MeshDecompositionSet(decomposition) + materialField.TypeSet(oc.FieldTypes.MATERIAL) + materialField.DecompositionSet(decomposition) materialField.GeometricFieldSet(geometricField) - materialField.VariableLabelSet(iron.FieldVariableTypes.U, "Material") + materialField.VariableLabelSet(oc.FieldVariableTypes.U, "Material") materialField.NumberOfVariablesSet(1) - materialField.NumberOfComponentsSet(iron.FieldVariableTypes.U,len(params)) - materialField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + materialField.NumberOfComponentsSet(oc.FieldVariableTypes.U,len(params)) + materialField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) if cellMLOption[0]: - print "----> CellML Material Field using gauss point interpolation <----\n" + print("----> CellML Material Field using gauss point interpolation <----\n") for component, param in enumerate(params, 1): - materialField.ComponentInterpolationSet(iron.FieldVariableTypes.U, component, - iron.FieldInterpolationTypes.GAUSS_POINT_BASED) + materialField.ComponentInterpolationSet(oc.FieldVariableTypes.U, component, + oc.FieldInterpolationTypes.GAUSS_POINT_BASED) else: - print "----> Material Field using constant interpolation <----\n" + print("----> Material Field using constant interpolation <----\n") for component, param in enumerate(params, 1): - materialField.ComponentInterpolationSet(iron.FieldVariableTypes.U, component, - iron.FieldInterpolationTypes.CONSTANT) + materialField.ComponentInterpolationSet(oc.FieldVariableTypes.U, component, + oc.FieldInterpolationTypes.CONSTANT) if option[0] == 2: # Tricubic Hermite if option[1] == 1: - materialField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + materialField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) elif option[1] == 2: - materialField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + materialField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) materialField.CreateFinish() for component, param in enumerate(params, 1): - materialField.ComponentValuesInitialiseDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, + materialField.ComponentValuesInitialiseDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, component, param) - materialField.ParameterSetUpdateStart(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES) - materialField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES) + materialField.ParameterSetUpdateStart(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES) + materialField.ParameterSetUpdateFinish(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES) # Output - print "----> Initialised " + str(len(params)) + " material parameters <----\n" + print("----> Initialised " + str(len(params)) + " material parameters <----\n") return materialField @@ -416,57 +440,57 @@ def MaterialFieldSetUp(materialFieldUserNumber, region, decomposition, geometric #=================================================================================# def DependentFieldSetUp(dependentFieldUserNumber, equationsSet, option, cellMLOption): # Set up dependent field - dependentField = iron.Field() + dependentField = oc.Field() equationsSet.DependentCreateStart(dependentFieldUserNumber, dependentField) - dependentField.VariableLabelSet(iron.FieldVariableTypes.U, "Dependent") + dependentField.VariableLabelSet(oc.FieldVariableTypes.U, "Dependent") if cellMLOption[0]: - print '----> Labelling dependent field strain and stress <----\n' - dependentField.VariableLabelSet(iron.FieldVariableTypes.U1, "Strain") - dependentField.VariableLabelSet(iron.FieldVariableTypes.U2, "Stress") + print('----> Labelling dependent field strain and stress <----\n') + dependentField.VariableLabelSet(oc.FieldVariableTypes.U1, "Strain") + dependentField.VariableLabelSet(oc.FieldVariableTypes.U2, "Stress") if option[0] == 1: # Trilinear for i in [1, 2, 3]: - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U, i, 1) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN, i, 1) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U, i, 1) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T, i, 1) - dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U, 4, - iron.FieldInterpolationTypes.ELEMENT_BASED) - dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.DELUDELN, 4, - iron.FieldInterpolationTypes.ELEMENT_BASED) + dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U, 4, + oc.FieldInterpolationTypes.ELEMENT_BASED) + dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.T, 4, + oc.FieldInterpolationTypes.ELEMENT_BASED) # Output - print "----> Use element based interpolation for hydrostatic pressure <----\n" + print("----> Use element based interpolation for hydrostatic pressure <----\n") elif option[0] == 2: # Tricubic Hermite for i in [1, 2, 3]: - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U, i, 1) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN, i, 1) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U, 4, 2) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN, 4, 2) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U, i, 1) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T, i, 1) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U, 4, 2) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T, 4, 2) - dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U, 4, - iron.FieldInterpolationTypes.NODE_BASED) - dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.DELUDELN, 4, - iron.FieldInterpolationTypes.NODE_BASED) + dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U, 4, + oc.FieldInterpolationTypes.NODE_BASED) + dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.T, 4, + oc.FieldInterpolationTypes.NODE_BASED) # Output - print "----> Interpolate hydrostatic pressure linearly <----\n" + print("----> Interpolate hydrostatic pressure linearly <----\n") if option[1] == 1: - dependentField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + dependentField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) # Output - print "----> Set up dependent field with unit scaling <----\n" + print("----> Set up dependent field with unit scaling <----\n") elif option[1] == 2: - dependentField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + dependentField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) # Output - print "----> Set up dependent field with arithmetic mean scaling <----\n" + print("----> Set up dependent field with arithmetic mean scaling <----\n") if cellMLOption[0]: for i in [1,2,3,4,5,6]: - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U1, i, 1) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U2, i, 1) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U1, i, 1) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U2, i, 1) equationsSet.DependentCreateFinish() @@ -481,22 +505,22 @@ def DependentFieldInitialise(dependentField, geometricField, hydroInit): # initial guess for hydrostatic pressure. # Copy over undeformed geometry to initialise dependent field. for i in [1, 2, 3]: - iron.Field.ParametersToFieldParametersComponentCopy(geometricField, iron.FieldVariableTypes.U, - iron.FieldParameterSetTypes.VALUES, i, dependentField, - iron.FieldVariableTypes.U, - iron.FieldParameterSetTypes.VALUES, i) + oc.Field.ParametersToFieldParametersComponentCopy(geometricField, oc.FieldVariableTypes.U, + oc.FieldParameterSetTypes.VALUES, i, dependentField, + oc.FieldVariableTypes.U, + oc.FieldParameterSetTypes.VALUES, i) # Output - print "----> Initialised dependent field with undeformed geometry <----\n" + print("----> Initialised dependent field with undeformed geometry <----\n") # Set hydrostatic pressure initial guess. - dependentField.ComponentValuesInitialise(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 4, + dependentField.ComponentValuesInitialise(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 4, hydroInit) - dependentField.ParameterSetUpdateStart(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES) - dependentField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES) + dependentField.ParameterSetUpdateStart(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES) + dependentField.ParameterSetUpdateFinish(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES) # Output - print "----> Initialised hydrostatic pressure guess of " + str(hydroInit) + " <----\n" + print("----> Initialised hydrostatic pressure guess of " + str(hydroInit) + " <----\n") #=================================================================================# @@ -511,38 +535,38 @@ def DependentFieldWarmStart(dependentField, deformedGeomDOFs, deformedHydro, opt for component in [1,2,3]: for node in range(1, numNodes+1): value = deformedGeomDOFs[component-1, node-1] - dependentField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, + dependentField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, 1, node, component, value) for e in range(1, len(deformedHydro)+1): - dependentField.ParameterSetUpdateElementDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, e, + dependentField.ParameterSetUpdateElementDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, e, 4, deformedHydro[e-1]) elif option[0] == 2: # Initialise dependent field to deformed warmstart solution. numNodes = len(deformedGeomDOFs[0,:,0]) for component in [1,2,3]: - print 'Component: ', component + print('Component: ', component) for node in range(1, numNodes+1): - print ' Node number: ', node + print(' Node number: ', node) for deriv in [1,2,3,4,5,6,7,8]: value = deformedGeomDOFs[component-1,node-1,deriv-1] - print ' value: ', value - dependentField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES, 1, + print(' value: ', value) + dependentField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES, 1, deriv, node, component, value) # Output - print "----> Initialised dependent field with warm-start geometry <----\n" + print("----> Initialised dependent field with warm-start geometry <----\n") # Set hydrostatic pressure initial guess. for node in range(1,numNodes+1): - dependentField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES,1, + dependentField.ParameterSetUpdateNodeDP(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES,1, 1, node, 4, deformedHydro[node-1]) - dependentField.ParameterSetUpdateStart(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES) - dependentField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES) + dependentField.ParameterSetUpdateStart(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES) + dependentField.ParameterSetUpdateFinish(oc.FieldVariableTypes.U, oc.FieldParameterSetTypes.VALUES) # Output - print "----> Initialised warm-start hydrostatic pressure of " + str(deformedHydro) + " <----\n" + print("----> Initialised warm-start hydrostatic pressure of " + str(deformedHydro) + " <----\n") #=================================================================================# @@ -560,7 +584,7 @@ def ParseWarmStart(filename, option): fid = open(filename, 'r') except IOError: - print 'ERROR: Unable to open ', filename + print('ERROR: Unable to open ', filename) return try: @@ -647,7 +671,7 @@ def ParseWarmStart(filename, option): temp = temp.split() num = temp[len(temp)-1] hydro = numpy.array(hydro) - print hydro + print(hydro) return nodes, hydro finally: fid.close() @@ -661,7 +685,7 @@ def ParseWarmStart(filename, option): fid = open(filename, 'r') except IOError: - print 'ERROR: Unable to open ', filename + print('ERROR: Unable to open ', filename) return try: @@ -742,7 +766,7 @@ def ParseWarmStart(filename, option): temp = temp.split() num = temp[len(temp)-1] hydro = numpy.array(hydro) - print hydro + print(hydro) return nodes, hydro finally: fid.close() @@ -758,13 +782,13 @@ def CellMLSetUp(cellMLUserNumber, cellMLModelsFieldUserNumber, cellMLParametersF # This function sets up the CellML environment for defining constitutive models. cellMLModelIndex = 1 - cellML = iron.CellML() + cellML = oc.CellML() cellML.CreateStart(cellMLUserNumber, region) cellML.ModelImport(filename) strain = ["E11", "E12", "E13", "E22", "E23", "E33"] stress2PK = ["Tdev11", "Tdev12", "Tdev13", "Tdev22", "Tdev23", "Tdev33"] - # Set strains as known in CellML. These will be fed into the model from iron. + # Set strains as known in CellML. These will be fed into the model from oc. for i in range(0, 6): cellML.VariableSetAsKnown(cellMLModelIndex, "equations/" + strain[i]) for component, parameter in enumerate(parameters): @@ -780,62 +804,62 @@ def CellMLSetUp(cellMLUserNumber, cellMLModelsFieldUserNumber, cellMLParametersF # Map the strain from dependentField U1 variable to CellML. for component, variable in enumerate(strain, 1): - #print "----> Mapping strain ", str(variable)+ " to CellML <----\n" - cellML.CreateFieldToCellMLMap(dependentField, iron.FieldVariableTypes.U1, component, - iron.FieldParameterSetTypes.VALUES, cellMLModelIndex, "equations/" + variable, - iron.FieldParameterSetTypes.VALUES) + #print("----> Mapping strain ", str(variable)+ " to CellML <----\n") + cellML.CreateFieldToCellMLMap(dependentField, oc.FieldVariableTypes.U1, component, + oc.FieldParameterSetTypes.VALUES, cellMLModelIndex, "equations/" + variable, + oc.FieldParameterSetTypes.VALUES) # Map the material parameters from material field to CellML. for component, parameter in enumerate(parameters, 1): - #print "----> Mapping parameter ", str(parameter)+ " to CellML <----\n" - cellML.CreateFieldToCellMLMap(materialField, iron.FieldVariableTypes.U, component, - iron.FieldParameterSetTypes.VALUES, cellMLModelIndex, "equations/" + parameter, - iron.FieldParameterSetTypes.VALUES) + #print("----> Mapping parameter ", str(parameter)+ " to CellML <----\n") + cellML.CreateFieldToCellMLMap(materialField, oc.FieldVariableTypes.U, component, + oc.FieldParameterSetTypes.VALUES, cellMLModelIndex, "equations/" + parameter, + oc.FieldParameterSetTypes.VALUES) # Map the stress from CellML to dependentFieldU2 variable for component, variable in enumerate(stress2PK, 1): - #print "----> Mapping stress ", str(variable)+ " to CellML <----\n" - cellML.CreateCellMLToFieldMap(cellMLModelIndex, "equations/" + variable, iron.FieldParameterSetTypes.VALUES, - dependentField, iron.FieldVariableTypes.U2, component, - iron.FieldParameterSetTypes.VALUES) + #print("----> Mapping stress ", str(variable)+ " to CellML <----\n") + cellML.CreateCellMLToFieldMap(cellMLModelIndex, "equations/" + variable, oc.FieldParameterSetTypes.VALUES, + dependentField, oc.FieldVariableTypes.U2, component, + oc.FieldParameterSetTypes.VALUES) cellML.FieldMapsCreateFinish() - print "----> Finished mapping variables to CellML <----\n" + print("----> Finished mapping variables to CellML <----\n") # Create models field for CellML - CellMLModelsField = iron.Field() + CellMLModelsField = oc.Field() cellML.ModelsFieldCreateStart(cellMLModelsFieldUserNumber, CellMLModelsField) if option[0] == 2: # Tricubic Hermite if option[1] == 1: - CellMLModelsField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + CellMLModelsField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) elif option[1] == 2: - CellMLModelsField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + CellMLModelsField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) cellML.ModelsFieldCreateFinish() - print "----> Finished creating models field for CellML <----\n" + print("----> Finished creating models field for CellML <----\n") # No need to create a state field since we aren't integrating. # Create parameters field for CellML, this is used as the strain field. - CellMLParametersField = iron.Field() + CellMLParametersField = oc.Field() cellML.ParametersFieldCreateStart(cellMLParametersFieldUserNumber, CellMLParametersField) if option[0] == 2: # Tricubic Hermite if option[1] == 1: - CellMLParametersField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + CellMLParametersField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) elif option[1] == 2: - CellMLParametersField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + CellMLParametersField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) cellML.ParametersFieldCreateFinish() - print "----> Finished creating parameters field for CellML <----\n" + print("----> Finished creating parameters field for CellML <----\n") # Create intermediate field for CellML, this is used as the stress field. - CellMLIntermediateField = iron.Field() + CellMLIntermediateField = oc.Field() cellML.IntermediateFieldCreateStart(cellMLIntermediateFieldUserNumber, CellMLIntermediateField) if option[0] == 2: # Tricubic Hermite if option[1] == 1: - CellMLIntermediateField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + CellMLIntermediateField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) elif option[1] == 2: - CellMLIntermediateField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + CellMLIntermediateField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) cellML.IntermediateFieldCreateFinish() - print "----> Finished creating intermediate field for CellML <----\n" + print("----> Finished creating intermediate field for CellML <----\n") return cellML, CellMLModelsField, CellMLParametersField, CellMLIntermediateField #=================================================================================# @@ -843,30 +867,30 @@ def CellMLSetUp(cellMLUserNumber, cellMLModelsFieldUserNumber, cellMLParametersF #=================================================================================# def StrainFieldSetUp(strainFieldUserNumber, region, decomposition, geometricField, equationsSet, option): # Set up strain field for output - strainField = iron.Field() + strainField = oc.Field() strainField.CreateStart(strainFieldUserNumber, region) - strainField.MeshDecompositionSet(decomposition) - strainField.TypeSet(iron.FieldTypes.GENERAL) + strainField.DecompositionSet(decomposition) + strainField.TypeSet(oc.FieldTypes.GENERAL) strainField.GeometricFieldSet(geometricField) - strainField.DependentTypeSet(iron.FieldDependentTypes.DEPENDENT) - strainField.VariableTypesSet([iron.FieldVariableTypes.U]) - strainField.VariableLabelSet(iron.FieldVariableTypes.U, "Strain") - strainField.NumberOfComponentsSet(iron.FieldVariableTypes.U, 6) + strainField.DependentTypeSet(oc.FieldDependentTypes.DEPENDENT) + strainField.VariableTypesSet([oc.FieldVariableTypes.U]) + strainField.VariableLabelSet(oc.FieldVariableTypes.U, "Strain") + strainField.NumberOfComponentsSet(oc.FieldVariableTypes.U, 6) for component in [1,2,3,4,5,6]: - strainField.ComponentInterpolationSet(iron.FieldVariableTypes.U, component, - iron.FieldInterpolationTypes.GAUSS_POINT_BASED) + strainField.ComponentInterpolationSet(oc.FieldVariableTypes.U, component, + oc.FieldInterpolationTypes.GAUSS_POINT_BASED) if option[0]==2: if option[1]==1: - strainField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + strainField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) elif option[1]==2: - strainField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + strainField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) strainField.CreateFinish() equationsSet.DerivedCreateStart(strainFieldUserNumber, strainField) - equationsSet.DerivedVariableSet(iron.EquationsSetDerivedTypes.STRAIN, iron.FieldVariableTypes.U) + equationsSet.DerivedVariableSet(oc.EquationsSetDerivedTypes.STRAIN, oc.FieldVariableTypes.U) equationsSet.DerivedCreateFinish() return strainField @@ -885,71 +909,71 @@ def matrixFromSymmetricComponents(components): def EquationsSetSetUp(equationsSet): # Set up standard options for problem and solvers. # Create equations - equations = iron.Equations() + equations = oc.Equations() equationsSet.EquationsCreateStart(equations) - equations.SparsityTypeSet(iron.EquationsSparsityTypes.SPARSE) - equations.OutputTypeSet(iron.EquationsOutputTypes.NONE) - #equations.OutputTypeSet(iron.EquationsOutputTypes.MATRIX) + equations.SparsityTypeSet(oc.EquationsSparsityTypes.SPARSE) + equations.OutputTypeSet(oc.EquationsOutputTypes.NONE) + #equations.OutputTypeSet(oc.EquationsOutputTypes.MATRIX) equationsSet.EquationsCreateFinish() #=================================================================================# #=================================================================================# -def ProblemSolverSetup(equationsSet,problemUserNumber,maxIter, TOL, cellMLOption): +def ProblemSolverSetup(context,equationsSet,problemUserNumber,maxIter, TOL, cellMLOption): # This function sets up the problem as well as the solver options. - print "----> Set up equations <----\n" + print("----> Set up equations <----\n") # Define problem - problem = iron.Problem() + problem = oc.Problem() if cellMLOption[0]: - problemSpecification = [iron.ProblemClasses.ELASTICITY, - iron.ProblemTypes.FINITE_ELASTICITY, - iron.ProblemSubtypes.FINITE_ELASTICITY_CELLML] + problemSpecification = [oc.ProblemClasses.ELASTICITY, + oc.ProblemTypes.FINITE_ELASTICITY, + oc.ProblemSubtypes.FINITE_ELASTICITY_CELLML] else: - problemSpecification = [iron.ProblemClasses.ELASTICITY, - iron.ProblemTypes.FINITE_ELASTICITY, - iron.ProblemSubtypes.NONE] + problemSpecification = [oc.ProblemClasses.ELASTICITY, + oc.ProblemTypes.FINITE_ELASTICITY, + oc.ProblemSubtypes.STATIC_FINITE_ELASTICITY] - problem.CreateStart(problemUserNumber, problemSpecification) + problem.CreateStart(problemUserNumber,context,problemSpecification) problem.CreateFinish() # Output - print "----> Set up problem <----\n" + print("----> Set up problem <----\n") # Create control loops problem.ControlLoopCreateStart() - controlLoop = iron.ControlLoop() - problem.ControlLoopGet([iron.ControlLoopIdentifiers.NODE], controlLoop) - #controlLoop.TypeSet(iron.ProblemControlLoopTypes.WHILE_LOOP) + controlLoop = oc.ControlLoop() + problem.ControlLoopGet([oc.ControlLoopIdentifiers.NODE], controlLoop) + controlLoop.TypeSet(oc.ControlLoopTypes.WHILE_LOOP) #controlLoop.IterationsSet(1,1,1) controlLoop.MaximumIterationsSet(maxIter) #controlLoop.MaximumIterationsSet(3) problem.ControlLoopCreateFinish() # Output - print "----> Set up control loop <----\n" + print("----> Set up control loop <----\n") # Create nonlinear numerical solver - linearSolver = iron.Solver() - nonLinearSolver = iron.Solver() + linearSolver = oc.Solver() + nonLinearSolver = oc.Solver() problem.SolversCreateStart() - problem.SolverGet([iron.ControlLoopIdentifiers.NODE], 1, nonLinearSolver) - nonLinearSolver.OutputTypeSet(iron.SolverOutputTypes.PROGRESS) - nonLinearSolver.NewtonJacobianCalculationTypeSet(iron.JacobianCalculationTypes.FD) + problem.SolverGet([oc.ControlLoopIdentifiers.NODE], 1, nonLinearSolver) + nonLinearSolver.OutputTypeSet(oc.SolverOutputTypes.PROGRESS) + nonLinearSolver.NewtonJacobianCalculationTypeSet(oc.JacobianCalculationTypes.FD) nonLinearSolver.NewtonAbsoluteToleranceSet(1e-12) nonLinearSolver.NewtonSolutionToleranceSet(1e-12) nonLinearSolver.NewtonRelativeToleranceSet(1e-12) - nonLinearSolver.NewtonConvergenceTestTypeSet(iron.NewtonConvergenceTypes.PETSC_DEFAULT) + nonLinearSolver.NewtonConvergenceTestTypeSet(oc.NewtonConvergenceTypes.PETSC_DEFAULT) nonLinearSolver.NewtonLinearSolverGet(linearSolver) - #nonLinearSolver.NewtonLineSearchTypeSet(iron.NewtonLineSearchTypes.LINEAR) + #nonLinearSolver.NewtonLineSearchTypeSet(oc.NewtonLineSearchTypes.LINEAR) #nonLinearSolver.NewtonLineSearchAlphaSet(1e-6) #nonLinearSolver.NewtonLineSearchMaxStepSet(1e5) #nonLinearSolver.NewtonLineSearchMonitorOutputSet() #nonLinearSolver.NewtonLineSearchStepTolSet(1e-5) - linearSolver.LinearTypeSet(iron.LinearSolverTypes.DIRECT) - #linearSolver.LinearDirectTypeSet(iron.DirectLinearSolverTypes.LU) - #linearSolver.LibraryTypeSet(iron.SolverLibraries.MUMPS) + linearSolver.LinearTypeSet(oc.LinearSolverTypes.DIRECT) + #linearSolver.LinearDirectTypeSet(oc.DirectLinearSolverTypes.LU) + #linearSolver.LibraryTypeSet(oc.SolverLibraries.MUMPS) problem.SolversCreateFinish() if cellMLOption[0]: - cellMLSolver = iron.Solver() - cellMLEquations = iron.CellMLEquations() + cellMLSolver = oc.Solver() + cellMLEquations = oc.CellMLEquations() problem.CellMLEquationsCreateStart() nonLinearSolver.NewtonCellMLSolverGet(cellMLSolver) cellMLSolver.CellMLEquationsGet(cellMLEquations) @@ -957,18 +981,18 @@ def ProblemSolverSetup(equationsSet,problemUserNumber,maxIter, TOL, cellMLOption problem.CellMLEquationsCreateFinish() # Output - print "----> Set up linear and nonlinear solvers <----\n" + print("----> Set up linear and nonlinear solvers <----\n") # Add solver equations sets which encompass the physics - solverEquations = iron.SolverEquations() - solver = iron.Solver() + solverEquations = oc.SolverEquations() + solver = oc.Solver() problem.SolverEquationsCreateStart() - problem.SolverGet([iron.ControlLoopIdentifiers.NODE], 1, solver) + problem.SolverGet([oc.ControlLoopIdentifiers.NODE], 1, solver) solver.SolverEquationsGet(solverEquations) - solverEquations.SparsityTypeSet(iron.SolverEquationsSparsityTypes.SPARSE) + solverEquations.SparsityTypeSet(oc.SolverEquationsSparsityTypes.SPARSE) equationSetIndex = solverEquations.EquationsSetAdd(equationsSet) problem.SolverEquationsCreateFinish() # Output - print "----> Set up solver with equations <----\n" + print("----> Set up solver with equations <----\n") return problem, solverEquations #=================================================================================# @@ -979,85 +1003,85 @@ def BCCubeSingleFace(solverEquations, dependentField, appliedFace, faceNormal, a # This function sets up the boundary conditions for dealing with BC's on a # single face of a cube. # Set up - boundaryConditions = iron.BoundaryConditions() + boundaryConditions = oc.BoundaryConditions() solverEquations.BoundaryConditionsCreateStart(boundaryConditions) # Initialise fixed faces node values. for node in fixXFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 1, - iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 1, + oc.BoundaryConditionsTypes.FIXED, 0.0) for node in fixYFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 2, - iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 2, + oc.BoundaryConditionsTypes.FIXED, 0.0) for node in fixZFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 3, - iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 3, + oc.BoundaryConditionsTypes.FIXED, 0.0) if option[0] == 2: # Fix derivatives if faceNormal == 1: - derivFix = [iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3] + derivFix = [oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S3] elif faceNormal == 2: - derivFix = [iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3] + derivFix = [oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S3] else: - derivFix = [iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2] + derivFix = [oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2] for node in range(1,numNodes+1): for j in derivFix: for component in [1,2,3]: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, j, node, - component, iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, j, node, + component, oc.BoundaryConditionsTypes.FIXED, 0.0) # Fix all second and third derivatives. for i in range(1, numNodes + 1): - for j in [iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3]: + for j in [oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3]: for k in [1, 2, 3]: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, j, i, k, - iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, j, i, k, + oc.BoundaryConditionsTypes.FIXED, 0.0) # Output - print "----> Implemented fixed boundary conditions <----\n" + print("----> Implemented fixed boundary conditions <----\n") # Initialise applied faces. if optionBC == 1: # Option 1: Compression/extension for node in appliedFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, - iron.BoundaryConditionsTypes.FIXED, increm) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, + oc.BoundaryConditionsTypes.FIXED, increm) # Output - print "----> Implemented compression/extension boundary condition of " + str(increm) + " <----\n" + print("----> Implemented compression/extension boundary condition of " + str(increm) + " <----\n") elif optionBC == 2: # Option 2: Force for node in appliedFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.DELUDELN, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, - iron.BoundaryConditionsTypes.FIXED_INCREMENTED, increm) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.T, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, + oc.BoundaryConditionsTypes.FIXED_INCREMENTED, increm) # Output - print "----> Implemented force boundary condition of " + str(increm) + "N <----\n" + print("----> Implemented force boundary condition of " + str(increm) + "N <----\n") elif optionBC == 3: # Option 3: Pressure for node in appliedFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.DELUDELN, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, faceNormal, - iron.BoundaryConditionsTypes.PRESSURE_INCREMENTED, increm) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.T, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, faceNormal, + oc.BoundaryConditionsTypes.PRESSURE_INCREMENTED, increm) # Output - print "----> Implemented pressure boundary condition of " + str(increm) + " kPa <----\n" + print("----> Implemented pressure boundary condition of " + str(increm) + " kPa <----\n") solverEquations.BoundaryConditionsCreateFinish() @@ -1069,79 +1093,79 @@ def BCCantilever(solverEquations, dependentField, appliedFace, faceNormal, appli fixBackFace, fixedFaceNormal, option): # This function sets up the BC for a cantilever problem. # Set up - boundaryConditions = iron.BoundaryConditions() + boundaryConditions = oc.BoundaryConditions() solverEquations.BoundaryConditionsCreateStart(boundaryConditions) # Initialise fixed faces node values. for component in [1, 2, 3]: for node in fixBackFace: for component in [1,2,3]: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, component, - iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, component, + oc.BoundaryConditionsTypes.FIXED, 0.0) if option[0] == 2: - # print 'Node number ', node + # print('Node number ', node) # Fix derivatives if fixedFaceNormal == 1: - #print "Fixed back normal is 1. " - derivFix = [iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3] + #print("Fixed back normal is 1. ") + derivFix = [oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3] elif fixedFaceNormal == 2: - derivFix = [iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3] + derivFix = [oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3] else: - derivFix = [iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3] + derivFix = [oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S2_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S2_S3] for deriv in derivFix: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, deriv, node, - component, iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, deriv, node, + component, oc.BoundaryConditionsTypes.FIXED, 0.0) # Output - print "----> Implemented fixed boundary conditions <----\n" + print("----> Implemented fixed boundary conditions <----\n") # Initialise applied faces. if optionBC == 1: # Option 1: Compression/extension for node in appliedFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, - iron.BoundaryConditionsTypes.FIXED, increm) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, + oc.BoundaryConditionsTypes.FIXED, increm) # Output - print "----> Implemented compression/extension boundary condition of " + str(increm) + " <----\n" + print("----> Implemented compression/extension boundary condition of " + str(increm) + " <----\n") elif optionBC == 2: # Option 2: Force for node in appliedFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.DELUDELN, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, - iron.BoundaryConditionsTypes.FIXED_INCREMENTED, increm) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.T, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, appliedDirection, + oc.BoundaryConditionsTypes.FIXED_INCREMENTED, increm) # Output - print "----> Implemented force boundary condition of " + str(increm) + "N <----\n" + print("----> Implemented force boundary condition of " + str(increm) + "N <----\n") elif optionBC == 3: # Option 3: Pressure - print 'Pressure applied on: ' + print('Pressure applied on: ') for node in appliedFace: - print 'Node ', node - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.DELUDELN, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, faceNormal, - iron.BoundaryConditionsTypes.PRESSURE_INCREMENTED, increm) - print 'Face normal ', faceNormal + print('Node ', node) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.T, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, faceNormal, + oc.BoundaryConditionsTypes.PRESSURE_INCREMENTED, increm) + print('Face normal ', faceNormal) # Output - print "----> Implemented pressure boundary condition of " + str(increm) + " kPa <----\n" + print("----> Implemented pressure boundary condition of " + str(increm) + " kPa <----\n") solverEquations.BoundaryConditionsCreateFinish() @@ -1153,37 +1177,37 @@ def BCEndoPressure(solverEquations, dependentField, endoFace, pressure, basalFac # This function sets up the BC for a LV inflation problem where pressure is applied # on the endocardial surface. # Set up - boundaryConditions = iron.BoundaryConditions() + boundaryConditions = oc.BoundaryConditions() solverEquations.BoundaryConditionsCreateStart(boundaryConditions) if option[0] == 1: - derivFix = [iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV] + derivFix = [oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV] else: - derivFix = [iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3, - iron.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3] + derivFix = [oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S3, + oc.GlobalDerivativeConstants.GLOBAL_DERIV_S1_S3] # Fix basal nodes and derivatives. for component in [1, 2, 3]: for node in basalFace: for deriv in derivFix: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, - deriv, node, component,iron.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, + deriv, node, component,oc.BoundaryConditionsTypes.FIXED, 0.0) # Apply pressure BC on endocardial nodes. for node in endoFace: - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.DELUDELN, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 3, - iron.BoundaryConditionsTypes.PRESSURE_INCREMENTED, pressure) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.T, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, node, 3, + oc.BoundaryConditionsTypes.PRESSURE_INCREMENTED, pressure) """ for component in [1,2,3]: - boundaryConditions.ConstrainNodeDofsEqual(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, component, apexEndoNodes) - boundaryConditions.ConstrainNodeDofsEqual(dependentField, iron.FieldVariableTypes.U, 1, - iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV, component, apexEpiNodes) + boundaryConditions.ConstrainNodeDofsEqual(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, component, apexEndoNodes) + boundaryConditions.ConstrainNodeDofsEqual(dependentField, oc.FieldVariableTypes.U, 1, + oc.GlobalDerivativeConstants.NO_GLOBAL_DERIV, component, apexEpiNodes) """ solverEquations.BoundaryConditionsCreateFinish() @@ -1194,32 +1218,32 @@ def BCEndoPressure(solverEquations, dependentField, endoFace, pressure, basalFac def ExportResults(dependentField, deformedFieldUserNumber, decomposition, region, filename, option): # This function exports the results of simulation to exnode and exelem files. # Copy over deformed field. - deformedField = iron.Field() + deformedField = oc.Field() deformedField.CreateStart(deformedFieldUserNumber, region) - deformedField.MeshDecompositionSet(decomposition) - deformedField.TypeSet(iron.FieldTypes.GEOMETRIC) - deformedField.VariableLabelSet(iron.FieldVariableTypes.U, "DeformedGeometry") + deformedField.DecompositionSet(decomposition) + deformedField.TypeSet(oc.FieldTypes.GEOMETRIC) + deformedField.VariableLabelSet(oc.FieldVariableTypes.U, "DeformedGeometry") if option[0] == 1: # Trilinear. for component in [1, 2, 3]: - deformedField.ComponentMeshComponentSet(iron.FieldVariableTypes.U, component, 1) - deformedField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + deformedField.ComponentMeshComponentSet(oc.FieldVariableTypes.U, component, 1) + deformedField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) elif option[0] == 2: # Tricubic hermite. Geometry interpolated using cubic hermite basis (2nd mesh component). for component in [1, 2, 3]: - deformedField.ComponentMeshComponentSet(iron.FieldVariableTypes.U, component, 1) + deformedField.ComponentMeshComponentSet(oc.FieldVariableTypes.U, component, 1) if option[1] == 1: - deformedField.ScalingTypeSet(iron.FieldScalingTypes.UNIT) + deformedField.ScalingTypeSet(oc.FieldScalingTypes.UNIT) elif option[1] == 2: - deformedField.ScalingTypeSet(iron.FieldScalingTypes.ARITHMETIC_MEAN) + deformedField.ScalingTypeSet(oc.FieldScalingTypes.ARITHMETIC_MEAN) deformedField.CreateFinish() for component in [1, 2, 3]: - dependentField.ParametersToFieldParametersComponentCopy(iron.FieldVariableTypes.U, - iron.FieldParameterSetTypes.VALUES, component, - deformedField, iron.FieldVariableTypes.U, - iron.FieldParameterSetTypes.VALUES, component) + dependentField.ParametersToFieldParametersComponentCopy(oc.FieldVariableTypes.U, + oc.FieldParameterSetTypes.VALUES, component, + deformedField, oc.FieldVariableTypes.U, + oc.FieldParameterSetTypes.VALUES, component) dependentField.Destroy() #deformedField.Destroy() # Export deformation. @@ -1227,14 +1251,14 @@ def ExportResults(dependentField, deformedFieldUserNumber, decomposition, region if not os.path.exists("./results"): os.makedirs("./results") - fields = iron.Fields() + fields = oc.Fields() fields.CreateRegion(region) fields.NodesExport("./results/" + filename, "FORTRAN") fields.ElementsExport("./results/" + filename, "FORTRAN") fields.Finalise() # Output - print "----> Export deformed geometric solutions <----\n" + print("----> Export deformed geometric solutions <----\n") #=================================================================================# @@ -1251,28 +1275,28 @@ def ExportStressStrain(elements, xiPositions, cellML, equationsSet, filename_dis file_disp = open(filename_disp, 'w') except IOError: - print 'ERROR: Unable to open ', filename_disp + print('ERROR: Unable to open ', filename_disp) return try: file_strain = open(filename_strain, 'w') except IOError: - print 'ERROR: Unable to open ', filename_strain + print('ERROR: Unable to open ', filename_strain) return try: file_stress2PK = open(filename_stress2PK, 'w') except IOError: - print 'ERROR: Unable to open ', filename_stress2PK + print('ERROR: Unable to open ', filename_stress2PK) return try: file_stressCauchy = open(filename_stressCauchy, 'w') except IOError: - print 'ERROR: Unable to open ', filename_stressCauchy + print('ERROR: Unable to open ', filename_stressCauchy) return # Write file headers for displacement @@ -1367,4 +1391,4 @@ def ExportStressStrain(elements, xiPositions, cellML, equationsSet, filename_dis file_stressCauchy.close() # Output - print "----> Export stresses and strains of deformed solution <----\n" + print("----> Export stresses and strains of deformed solution <----\n")