/FiniteElasticity/Cantilever/src/CantileverExample.py
Python | 346 lines | 221 code | 41 blank | 84 comment | 14 complexity | ce79b12173645ceea556e8f413d47c5e MD5 | raw file
1#> \file 2#> \author Chris Bradley 3#> \brief This is an example script to solve a finite elasticity equation using openCMISS calls in python. 4#> 5#> \section LICENSE 6#> 7#> Version: MPL 1.1/GPL 2.0/LGPL 2.1 8#> 9#> The contents of this file are subject to the Mozilla Public License 10#> Version 1.1 (the "License"); you may not use this file except in 11#> compliance with the License. You may obtain a copy of the License at 12#> http://www.mozilla.org/MPL/ 13#> 14#> Software distributed under the License is distributed on an "AS IS" 15#> basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the 16#> License for the specific language governing rights and limitations 17#> under the License. 18#> 19#> The Original Code is openCMISS 20#> 21#> The Initial Developer of the Original Code is University of Auckland, 22#> Auckland, New Zealand and University of Oxford, Oxford, United 23#> Kingdom. Portions created by the University of Auckland and University 24#> of Oxford are Copyright (C) 2007 by the University of Auckland and 25#> the University of Oxford. All Rights Reserved. 26#> 27#> Contributor(s): 28#> 29#> Alternatively, the contents of this file may be used under the terms of 30#> either the GNU General Public License Version 2 or later (the "GPL"), or 31#> the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 32#> in which case the provisions of the GPL or the LGPL are applicable instead 33#> of those above. if you wish to allow use of your version of this file only 34#> under the terms of either the GPL or the LGPL, and not to allow others to 35#> use your version of this file under the terms of the MPL, indicate your 36#> decision by deleting the provisions above and replace them with the notice 37#> and other provisions required by the GPL or the LGPL. if you do not delete 38#> the provisions above, a recipient may use your version of this file under 39#> the terms of any one of the MPL, the GPL or the LGPL. 40#> 41 42#> \example FiniteElasticity/Cantilever/src/CantileverExample.py 43## Example script to solve a finite elasticity equation using openCMISS calls in python. 44## \par Latest Builds: 45## \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/FiniteElasticity/Cantilever/build-intel'>Linux Intel Build</a> 46## \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/FiniteElasticity/Cantilever/build-gnu'>Linux GNU Build</a> 47#< 48 49#> Main script 50# Add Python bindings directory to PATH 51import sys, os 52 53sys.path.append(os.sep.join((os.environ['OPENCMISS_ROOT'],'cm','bindings','python'))) 54 55# Intialise OpenCMISS 56from opencmiss import CMISS 57 58# Set problem parameters 59 60width = 60.0 61length = 40.0 62height = 40.0 63density=9.0E-4 #in g mm^-3 64gravity=[0.0,0.0,-9.81] #in m s^-2 65 66UsePressureBasis = False 67NumberOfGaussXi = 2 68 69coordinateSystemUserNumber = 1 70regionUserNumber = 1 71basisUserNumber = 1 72pressureBasisUserNumber = 2 73generatedMeshUserNumber = 1 74meshUserNumber = 1 75decompositionUserNumber = 1 76geometricFieldUserNumber = 1 77fibreFieldUserNumber = 2 78materialFieldUserNumber = 3 79dependentFieldUserNumber = 4 80sourceFieldUserNumber = 5 81equationsSetFieldUserNumber = 6 82equationsSetUserNumber = 1 83problemUserNumber = 1 84 85CMISS.ErrorHandlingModeSet(CMISS.ErrorHandlingModes.TrapError) 86 87# Set all diganostic levels on for testing 88#CMISS.DiagnosticsSetOn(CMISS.DiagnosticTypes.All,[1,2,3,4,5],"Diagnostics",["DOMAIN_MAPPINGS_LOCAL_FROM_GLOBAL_CALCULATE"]) 89 90numberOfLoadIncrements = 2 91numberGlobalXElements = 1 92numberGlobalYElements = 1 93numberGlobalZElements = 1 94InterpolationType = 1 95if(numberGlobalZElements==0): 96 numberOfXi = 2 97else: 98 numberOfXi = 3 99 100# Get the number of computational nodes and this computational node number 101numberOfComputationalNodes = CMISS.ComputationalNumberOfNodesGet() 102computationalNodeNumber = CMISS.ComputationalNodeNumberGet() 103 104# Create a 3D rectangular cartesian coordinate system 105coordinateSystem = CMISS.CoordinateSystem() 106coordinateSystem.CreateStart(coordinateSystemUserNumber) 107coordinateSystem.DimensionSet(3) 108coordinateSystem.CreateFinish() 109 110# Create a region and assign the coordinate system to the region 111region = CMISS.Region() 112region.CreateStart(regionUserNumber,CMISS.WorldRegion) 113#region.LabelSet("Region") 114region.coordinateSystem = coordinateSystem 115region.CreateFinish() 116 117# Define basis 118basis = CMISS.Basis() 119basis.CreateStart(basisUserNumber) 120if InterpolationType == (1,2,3,4): 121 basis.type = CMISS.BasisTypes.LagrangeHermiteTP 122elif InterpolationType == (7,8,9): 123 basis.type = CMISS.BasisTypes.BasisSimplexType 124basis.numberOfXi = numberOfXi 125basis.interpolationXi = [CMISS.BasisInterpolationSpecifications.LinearLagrange]*numberOfXi 126if(NumberOfGaussXi>0): 127 basis.quadratureNumberOfGaussXi = [NumberOfGaussXi]*numberOfXi 128basis.CreateFinish() 129 130if(UsePressureBasis): 131 # Define pressure basis 132 pressureBasis = CMISS.Basis() 133 pressureBasis.CreateStart(pressureBasisUserNumber) 134 if InterpolationType == (1,2,3,4): 135 pressureBasis.type = CMISS.BasisTypes.LagrangeHermiteTP 136 elif InterpolationType == (7,8,9): 137 pressureBasis.type = CMISS.BasisTypes.BasisSimplexType 138 pressureBasis.numberOfXi = numberOfXi 139 pressureBasis.interpolationXi = [CMISS.BasisInterpolationSpecifications.LinearLagrange]*numberOfXi 140 if(NumberOfGaussXi>0): 141 pressureBasis.quadratureNumberOfGaussXi = [NumberOfGaussXi]*numberOfXi 142 pressureBasis.CreateFinish() 143 144# Start the creation of a generated mesh in the region 145generatedMesh = CMISS.GeneratedMesh() 146generatedMesh.CreateStart(generatedMeshUserNumber,region) 147generatedMesh.type = CMISS.GeneratedMeshTypes.Regular 148if(UsePressureBasis): 149 generatedMesh.basis = [basis,pressureBasis] 150else: 151 generatedMesh.basis = [basis] 152if(numberGlobalZElements==0): 153 generatedMesh.extent = [width,height] 154 generatedMesh.numberOfElements = [numberGlobalXElements,numberGlobalYElements] 155else: 156 generatedMesh.extent = [width,length,height] 157 generatedMesh.numberOfElements = [numberGlobalXElements,numberGlobalYElements,numberGlobalZElements] 158# Finish the creation of a generated mesh in the region 159mesh = CMISS.Mesh() 160generatedMesh.CreateFinish(meshUserNumber,mesh) 161 162# Create a decomposition for the mesh 163decomposition = CMISS.Decomposition() 164decomposition.CreateStart(decompositionUserNumber,mesh) 165decomposition.type = CMISS.DecompositionTypes.Calculated 166decomposition.numberOfDomains = numberOfComputationalNodes 167decomposition.CreateFinish() 168 169# Create a field for the geometry 170geometricField = CMISS.Field() 171geometricField.CreateStart(geometricFieldUserNumber,region) 172geometricField.MeshDecompositionSet(decomposition) 173geometricField.TypeSet(CMISS.FieldTypes.Geometric) 174geometricField.VariableLabelSet(CMISS.FieldVariableTypes.U,"Geometry") 175geometricField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.U,1,1) 176geometricField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.U,2,1) 177geometricField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.U,3,1) 178if InterpolationType == 4: 179 geometricField.fieldScalingType = CMISS.FieldScalingTypes.ArithmeticMean 180geometricField.CreateFinish() 181 182# Update the geometric field parameters from generated mesh 183CMISS.GeneratedMeshGeometricParametersCalculate(geometricField,generatedMesh) 184 185# Create a fibre field and attach it to the geometric field 186fibreField = CMISS.Field() 187fibreField.CreateStart(fibreFieldUserNumber,region) 188fibreField.TypeSet(CMISS.FieldTypes.Fibre) 189fibreField.MeshDecompositionSet(decomposition) 190fibreField.GeometricFieldSet(geometricField) 191fibreField.VariableLabelSet(CMISS.FieldVariableTypes.U,"Fibre") 192if InterpolationType == 4: 193 fibreField.fieldScalingType = CMISS.FieldScalingTypes.ArithmeticMean 194fibreField.CreateFinish() 195 196# Create the equations_set 197equationsSetField = CMISS.Field() 198equationsSet = CMISS.EquationsSet() 199equationsSet.CreateStart(equationsSetUserNumber,region,fibreField, \ 200 CMISS.EquationsSetClasses.Elasticity, 201 CMISS.EquationsSetTypes.FiniteElasticity, \ 202 CMISS.EquationsSetSubtypes.MooneyRivlin, \ 203 equationsSetFieldUserNumber, equationsSetField) 204equationsSet.CreateFinish() 205 206# Create the dependent field 207dependentField = CMISS.Field() 208equationsSet.DependentCreateStart(dependentFieldUserNumber,dependentField) 209dependentField.VariableLabelSet(CMISS.FieldVariableTypes.U,"Dependent") 210dependentField.ComponentInterpolationSet(CMISS.FieldVariableTypes.U,4,CMISS.FieldInterpolationTypes.ElementBased) 211dependentField.ComponentInterpolationSet(CMISS.FieldVariableTypes.DelUDelN,4,CMISS.FieldInterpolationTypes.ElementBased) 212if(UsePressureBasis): 213 # Set the pressure to be nodally based and use the second mesh component 214 if InterpolationType == 4: 215 dependentField.ComponentInterpolationSet(CMISS.FieldVariableTypes.U,4,CMISS.FieldInterpolationTypes.NodeBased) 216 dependentField.ComponentInterpolationSet(CMISS.FieldVariableTypes.DelUDelN,4,CMISS.FieldInterpolationTypes.NodeBased) 217 dependentField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.U,4,2) 218 dependentField.ComponentMeshComponentSet(CMISS.FieldVariableTypes.DelUDelN,4,2) 219if InterpolationType == 4: 220 dependentField.fieldScalingType = CMISS.FieldScalingTypes.ArithmeticMean 221equationsSet.DependentCreateFinish() 222 223 224# Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure 225CMISS.Field.ParametersToFieldParametersComponentCopy( \ 226 geometricField,CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,1, \ 227 dependentField,CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,1) 228CMISS.Field.ParametersToFieldParametersComponentCopy( \ 229 geometricField,CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,2, \ 230 dependentField,CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,2) 231CMISS.Field.ParametersToFieldParametersComponentCopy( \ 232 geometricField,CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,3, \ 233 dependentField,CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,3) 234CMISS.Field.ComponentValuesInitialiseDP( \ 235 dependentField,CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,4,-8.0) 236 237# Create the material field 238materialField = CMISS.Field() 239equationsSet.MaterialsCreateStart(materialFieldUserNumber,materialField) 240materialField.VariableLabelSet(CMISS.FieldVariableTypes.U,"Material") 241materialField.VariableLabelSet(CMISS.FieldVariableTypes.V,"Density") 242equationsSet.MaterialsCreateFinish() 243 244# Set Mooney-Rivlin constants c10 and c01 respectively. 245materialField.ComponentValuesInitialiseDP( \ 246 CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,1,2.0) 247materialField.ComponentValuesInitialiseDP( \ 248 CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,2,-14.0) 249materialField.ComponentValuesInitialiseDP( \ 250 CMISS.FieldVariableTypes.V,CMISS.FieldParameterSetTypes.FieldValues,1,density) 251 252#Create the source field with the gravity vector 253sourceField = CMISS.Field() 254equationsSet.SourceCreateStart(sourceFieldUserNumber,sourceField) 255sourceField.fieldScalingType = CMISS.FieldScalingTypes.ArithmeticMean 256if InterpolationType == 4: 257 sourceField.fieldScalingType = CMISS.FieldScalingTypes.ArithmeticMean 258else: 259 sourceField.fieldScalingType = CMISS.FieldScalingTypes.Unit 260equationsSet.SourceCreateFinish() 261 262#Set the gravity vector component values 263sourceField.ComponentValuesInitialiseDP( \ 264 CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,1,gravity[0]) 265sourceField.ComponentValuesInitialiseDP( \ 266 CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,2,gravity[1]) 267sourceField.ComponentValuesInitialiseDP( \ 268 CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.FieldValues,3,gravity[2]) 269 270# Create equations 271equations = CMISS.Equations() 272equationsSet.EquationsCreateStart(equations) 273equations.sparsityType = CMISS.EquationsSparsityTypes.Sparse 274equations.outputType = CMISS.EquationsOutputTypes.NONE 275equationsSet.EquationsCreateFinish() 276 277# Define the problem 278problem = CMISS.Problem() 279problem.CreateStart(problemUserNumber) 280problem.SpecificationSet(CMISS.ProblemClasses.Elasticity, \ 281 CMISS.ProblemTypes.FiniteElasticity, \ 282 CMISS.ProblemSubTypes.NONE) 283problem.CreateFinish() 284 285# Create the problem control loop 286problem.ControlLoopCreateStart() 287controlLoop = CMISS.ControlLoop() 288problem.ControlLoopGet([CMISS.ControlLoopIdentifiers.Node],controlLoop) 289controlLoop.MaximumIterationsSet(numberOfLoadIncrements) 290problem.ControlLoopCreateFinish() 291 292# Create problem solver 293nonLinearSolver = CMISS.Solver() 294linearSolver = CMISS.Solver() 295problem.SolversCreateStart() 296problem.SolverGet([CMISS.ControlLoopIdentifiers.Node],1,nonLinearSolver) 297nonLinearSolver.outputType = CMISS.SolverOutputTypes.Progress 298nonLinearSolver.NewtonJacobianCalculationTypeSet(CMISS.JacobianCalculationTypes.AnalyticCalculated) 299nonLinearSolver.NewtonLinearSolverGet(linearSolver) 300linearSolver.linearType = CMISS.LinearSolverTypes.Direct 301#linearSolver.libraryType = CMISS.SolverLibraries.LAPACK 302problem.SolversCreateFinish() 303 304# Create solver equations and add equations set to solver equations 305solver = CMISS.Solver() 306solverEquations = CMISS.SolverEquations() 307problem.SolverEquationsCreateStart() 308problem.SolverGet([CMISS.ControlLoopIdentifiers.Node],1,solver) 309solver.SolverEquationsGet(solverEquations) 310solverEquations.sparsityType = CMISS.SolverEquationsSparsityTypes.Sparse 311equationsSetIndex = solverEquations.EquationsSetAdd(equationsSet) 312problem.SolverEquationsCreateFinish() 313 314# Prescribe boundary conditions (absolute nodal parameters) 315boundaryConditions = CMISS.BoundaryConditions() 316solverEquations.BoundaryConditionsCreateStart(boundaryConditions) 317# Set x=0 nodes to no x displacment 318boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,1,1,CMISS.BoundaryConditionsTypes.Fixed,0.0) 319boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,3,1,CMISS.BoundaryConditionsTypes.Fixed,0.0) 320boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,5,1,CMISS.BoundaryConditionsTypes.Fixed,0.0) 321boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,7,1,CMISS.BoundaryConditionsTypes.Fixed,0.0) 322 323# Set y=0 nodes to no y displacement 324boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,1,2,CMISS.BoundaryConditionsTypes.Fixed,0.0) 325boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,3,2,CMISS.BoundaryConditionsTypes.Fixed,0.0) 326boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,5,2,CMISS.BoundaryConditionsTypes.Fixed,0.0) 327boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,7,2,CMISS.BoundaryConditionsTypes.Fixed,0.0) 328 329# Set z=0 nodes to no y displacement 330boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,1,3,CMISS.BoundaryConditionsTypes.Fixed,0.0) 331boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,3,3,CMISS.BoundaryConditionsTypes.Fixed,0.0) 332boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,5,3,CMISS.BoundaryConditionsTypes.Fixed,0.0) 333boundaryConditions.AddNode(dependentField,CMISS.FieldVariableTypes.U,1,1,7,3,CMISS.BoundaryConditionsTypes.Fixed,0.0) 334solverEquations.BoundaryConditionsCreateFinish() 335 336# Solve the problem 337problem.Solve() 338 339# Export results 340fields = CMISS.Fields() 341CMISS.FieldsTypeCreateRegion(region,fields) 342CMISS.FieldIONodesExport(fields,"../Cantilever","FORTRAN") 343CMISS.FieldIOElementsExport(fields,"../Cantilever","FORTRAN") 344fields.Finalise() 345 346