PageRenderTime 117ms CodeModel.GetById 16ms app.highlight 43ms RepoModel.GetById 1ms app.codeStats 0ms

/FiniteElasticity/Cantilever/src/CantileverExample.py

http://github.com/PrasadBabarendaGamage/examples
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