PageRenderTime 107ms CodeModel.GetById 24ms app.highlight 77ms RepoModel.GetById 1ms app.codeStats 0ms

/FiniteElasticity/QuadraticEllipsoid/src/QuadraticEllipsoidExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 559 lines | 382 code | 71 blank | 106 comment | 0 complexity | 85f98210a478311aca59b9681ddd0f59 MD5 | raw file
  1!> \file
  2!> \author Chris Bradley
  3!> \brief This is an example program to solve a finite elasticity equation using openCMISS calls.
  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): Jack Lee
 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/UniAxialExtension/src/UniAxialExtensionExample.f90
 43!! Example program to solve a finite elasticity equation using openCMISS calls.
 44!! \par Latest Builds:
 45!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/FiniteElasticity/UniAxialExtension/build-intel'>Linux Intel Build</a>
 46!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/FiniteElasticity/UniAxialExtension/build-gnu'>Linux GNU Build</a>
 47!<
 48
 49!> Main program
 50PROGRAM QUADRATICELLIPSOIDEEXAMPLE
 51
 52  USE OPENCMISS
 53  USE MPI
 54
 55#ifdef WIN32
 56  USE IFQWIN
 57#endif
 58
 59  IMPLICIT NONE
 60
 61  !Test program parameters
 62  
 63  REAL(CMISSDP), PARAMETER :: PI=3.14159_CMISSDP
 64  REAL(CMISSDP), PARAMETER :: LONG_AXIS=2.0_CMISSDP
 65  REAL(CMISSDP), PARAMETER :: SHORT_AXIS=1.0_CMISSDP
 66  REAL(CMISSDP), PARAMETER :: WALL_THICKNESS=0.5_CMISSDP
 67  REAL(CMISSDP), PARAMETER :: CUTOFF_ANGLE=1.5708_CMISSDP
 68  REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_INTERSECTION=1.73205_CMISSDP !Slope of fibres in base endocardium = 60 degrees
 69  REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_CHANGE=-3.4641_CMISSDP !Slope change of fibres from 60 to -60 degrees in transmural direction 
 70  REAL(CMISSDP), PARAMETER :: SHEET_SLOPE_BASE_ENDO=1.0_CMISSDP !Slope of sheet at base endocardium 
 71
 72  REAL(CMISSDP), PARAMETER :: INNER_PRESSURE=2.0_CMISSDP  !Positive is compressive
 73  REAL(CMISSDP), PARAMETER :: OUTER_PRESSURE=0.0_CMISSDP  !Positive is compressive
 74
 75  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalXElements=4  ! X ==NUMBER_GLOBAL_CIRCUMFERENTIAL_ELEMENTS
 76  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalYElements=4  ! Y ==NUMBER_GLOBAL_LONGITUDINAL_ELEMENTS
 77  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalZElements=1  ! Z ==NUMBER_GLOBAL_TRANSMURAL_ELEMENTS
 78  INTEGER(CMISSIntg) :: NumberOfDomains
 79
 80  INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
 81  INTEGER(CMISSIntg), PARAMETER :: NumberOfSpatialCoordinates=3
 82  INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=1
 83  INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=1
 84  INTEGER(CMISSIntg), PARAMETER :: QuadraticCollapsedBasisUserNumber=2
 85  INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=3
 86  INTEGER(CMISSIntg), PARAMETER :: LinearCollapsedBasisUserNumber=4
 87  INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=1
 88  INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=2
 89  INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=1
 90  INTEGER(CMISSIntg), PARAMETER :: DerivativeUserNumber=1
 91  
 92  INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshDimensions=3
 93  INTEGER(CMISSIntg), PARAMETER :: NumberOfXiCoordinates=3
 94  INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshComponents=2
 95  INTEGER(CMISSIntg), PARAMETER :: QuadraticMeshComponentNumber=1
 96  INTEGER(CMISSIntg), PARAMETER :: LinearMeshComponentNumber=2
 97  INTEGER(CMISSIntg), PARAMETER :: TotalNumberOfElements=1
 98
 99  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumber=1
100  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1
101  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
102
103  INTEGER(CMISSIntg), PARAMETER :: FieldFibreUserNumber=2
104  INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfVariables=1
105  INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfComponents=3
106
107  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialUserNumber=3
108  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfVariables=1
109  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfComponents=2
110
111  INTEGER(CMISSIntg), PARAMETER :: FieldDependentUserNumber=4
112  INTEGER(CMISSIntg), PARAMETER :: FieldDependentNumberOfVariables=2
113  INTEGER(CMISSIntg), PARAMETER :: FieldDependentNumberOfComponents=4
114
115  INTEGER(CMISSIntg), PARAMETER :: EquationSetUserNumber=1
116  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumber=13
117  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=1
118
119  !Program types
120
121
122  !Program variables
123
124  INTEGER(CMISSIntg) :: MPI_IERROR
125  INTEGER(CMISSIntg) :: EquationsSetIndex  
126  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber
127  REAL(CMISSDP) :: FibreFieldAngle(3) 
128  REAL(CMISSDP) :: nu,theta,omega,XI3,XI3delta,XI2delta, zero
129  INTEGER(CMISSIntg) ::i,j,k,component_idx,node_idx,TOTAL_NUMBER_NODES_XI(3)
130  !For grabbing surfaces
131  INTEGER(CMISSIntg) :: InnerNormalXi,OuterNormalXi,TopNormalXi
132  INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodes(:)
133  INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodes(:)
134  INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodes(:)
135  INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
136  REAL(CMISSDP) :: XCoord,YCoord,ZCoord
137  LOGICAL :: X_FIXED,Y_FIXED,X_OKAY,Y_OKAY
138
139  !CMISS variables
140
141  TYPE(CMISSBasisType) :: QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis
142  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditions
143  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem, WorldCoordinateSystem
144  TYPE(CMISSMeshType) :: Mesh
145  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
146  TYPE(CMISSDecompositionType) :: Decomposition
147  TYPE(CMISSEquationsType) :: Equations
148  TYPE(CMISSEquationsSetType) :: EquationsSet
149  TYPE(CMISSFieldType) :: GeometricField,FibreField,MaterialField,DependentField,EquationsSetField
150  TYPE(CMISSFieldsType) :: Fields
151  TYPE(CMISSProblemType) :: Problem
152  TYPE(CMISSRegionType) :: Region,WorldRegion
153  TYPE(CMISSSolverType) :: Solver,LinearSolver
154  TYPE(CMISSSolverEquationsType) :: SolverEquations
155
156#ifdef WIN32
157  !Quickwin type
158  LOGICAL :: QUICKWIN_STATUS=.FALSE.
159  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
160#endif
161
162  !Generic CMISS variables
163  INTEGER(CMISSIntg) :: Err
164
165#ifdef WIN32
166  !Initialise QuickWin
167  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
168  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
169  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
170  !Set the window parameters
171  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
172  !If attempt fails set with system estimated values
173  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
174#endif
175
176  !Intialise cmiss
177  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
178
179  CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
180
181  WRITE(*,'(A)') "Program starting."
182
183  !Set all diganostic levels on for testing
184  CALL CMISSDiagnosticsSetOn(CMISS_FROM_DIAG_TYPE,(/1,2,3,4,5/),"Diagnostics",(/"PROBLEM_FINITE_ELEMENT_CALCULATE"/),Err)
185
186  !Get the number of computational nodes and this computational node number
187  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
188  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
189
190  !Broadcast the number of elements in the X,Y and Z directions and the number of partitions to the other computational nodes
191!   CALL MPI_BCAST(NumberGlobalXElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
192!   CALL MPI_BCAST(NumberGlobalYElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
193!   CALL MPI_BCAST(NumberGlobalZElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
194!   CALL MPI_BCAST(NumberOfDomains,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
195  NumberOfDomains=NumberOfComputationalNodes
196
197  !Create a CS - default is 3D rectangular cartesian CS with 0,0,0 as origin
198  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
199  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
200  CALL CMISSCoordinateSystem_TypeSet(CoordinateSystem,CMISS_COORDINATE_RECTANGULAR_CARTESIAN_TYPE,Err)
201  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NumberOfSpatialCoordinates,Err)
202  CALL CMISSCoordinateSystem_OriginSet(CoordinateSystem,(/0.0_CMISSDP,0.0_CMISSDP,0.0_CMISSDP/),Err)
203  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
204
205  !Create a region and assign the CS to the region
206  CALL CMISSRegion_Initialise(Region,Err)
207  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
208  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
209  CALL CMISSRegion_CreateFinish(Region,Err)
210
211  !Define basis functions - tri-linear Lagrange and tri-Quadratic Lagrange, each with collapsed variant
212    !Quadratic Basis
213  CALL CMISSBasis_Initialise(QuadraticBasis,Err)
214  CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
215  CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
216    & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
217  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
218    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
219  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err) !Have to do this
220  CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
221
222    !Collapsed Quadratic Basis
223  CALL CMISSBasis_Initialise(QuadraticCollapsedBasis,Err)
224  CALL CMISSBasis_CreateStart(QuadraticCollapsedBasisUserNumber,QuadraticCollapsedBasis,Err)
225  CALL CMISSBasis_TypeSet(QuadraticCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
226  CALL CMISSBasis_NumberOfXiSet(QuadraticCollapsedBasis,NumberOfXiCoordinates,Err)
227  CALL CMISSBasis_InterpolationXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
228       & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
229  CALL CMISSBasis_CollapsedXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED, &
230       & CMISS_BASIS_COLLAPSED_AT_XI0,CMISS_BASIS_NOT_COLLAPSED/),Err)
231  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticCollapsedBasis, &
232       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)  
233  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticCollapsedBasis,.true.,Err) !Have to do this
234  CALL CMISSBasis_CreateFinish(QuadraticCollapsedBasis,Err)
235
236    !Linear Basis
237  CALL CMISSBasis_Initialise(LinearBasis,Err)
238  CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
239  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
240    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
241  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
242  CALL CMISSBasis_CreateFinish(LinearBasis,Err)
243
244    !Collapsed Linear Basis
245  CALL CMISSBasis_Initialise(LinearCollapsedBasis,Err)
246  CALL CMISSBasis_CreateStart(LinearCollapsedBasisUserNumber,LinearCollapsedBasis,Err)
247  CALL CMISSBasis_TypeSet(LinearCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
248  CALL CMISSBasis_NumberOfXiSet(LinearCollapsedBasis,NumberOfXiCoordinates,Err)
249  CALL CMISSBasis_InterpolationXiSet(LinearCollapsedBasis,(/CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION, &
250       & CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION,CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION/),Err)
251  CALL CMISSBasis_CollapsedXiSet(LinearCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED,CMISS_BASIS_COLLAPSED_AT_XI0, &
252    & CMISS_BASIS_NOT_COLLAPSED/),Err)
253  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearCollapsedBasis, &
254       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
255  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearCollapsedBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
256  CALL CMISSBasis_CreateFinish(LinearCollapsedBasis,Err)
257
258  !Start the creation of a generated ellipsoid mesh
259  CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
260  CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
261  !Set up an ellipsoid mesh
262  CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_MESH_TYPE,Err)
263  !Set the quadratic and linear bases
264  CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,[QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis],Err)
265  !Define the mesh on the region
266  CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/LONG_AXIS,SHORT_AXIS,WALL_THICKNESS,CUTOFF_ANGLE/),Err)
267  CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NumberGlobalXElements,NumberGlobalYElements, &
268    & NumberGlobalZElements/),Err)
269  
270  !Finish the creation of a generated mesh in the region
271  CALL CMISSMesh_Initialise(Mesh,Err)
272  CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
273
274  !Create a decomposition
275  CALL CMISSDecomposition_Initialise(Decomposition,Err)
276  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
277  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
278  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
279  CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.TRUE.,Err)
280  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
281
282  !Create a field to put the geometry (default is geometry)
283  CALL CMISSField_Initialise(GeometricField,Err)
284  CALL CMISSField_CreateStart(FieldGeometryUserNumber,Region,GeometricField,Err)
285  CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
286  CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)  
287  CALL CMISSField_NumberOfVariablesSet(GeometricField,FieldGeometryNumberOfVariables,Err)
288  CALL CMISSField_NumberOfComponentsSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)  
289  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
290  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
291  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
292  CALL CMISSField_CreateFinish(GeometricField,Err)
293
294  !Update the geometric field parameters
295  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricField,Err)
296
297  !Create a fibre field and attach it to the geometric field  
298  CALL CMISSField_Initialise(FibreField,Err)
299  CALL CMISSField_CreateStart(FieldFibreUserNumber,Region,FibreField,Err)
300  CALL CMISSField_TypeSet(FibreField,CMISS_FIELD_FIBRE_TYPE,Err)
301  CALL CMISSField_MeshDecompositionSet(FibreField,Decomposition,Err)        
302  CALL CMISSField_GeometricFieldSet(FibreField,GeometricField,Err)
303  CALL CMISSField_NumberOfVariablesSet(FibreField,FieldFibreNumberOfVariables,Err)
304  CALL CMISSField_NumberOfComponentsSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreNumberOfComponents,Err)  
305  CALL CMISSField_ComponentMeshComponentSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
306  CALL CMISSField_ComponentMeshComponentSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
307  CALL CMISSField_ComponentMeshComponentSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
308  CALL CMISSField_CreateFinish(FibreField,Err)
309
310  !Set Fibre directions (this block is parallel-untested)
311  node_idx=0  
312  !This is valid only for quadratic basis functions
313  TOTAL_NUMBER_NODES_XI(1)=NumberGlobalXElements*2
314  TOTAL_NUMBER_NODES_XI(2)=NumberGlobalYElements*2+1
315  TOTAL_NUMBER_NODES_XI(3)=NumberGlobalZElements*2+1
316   
317  XI2delta=(PI-CUTOFF_ANGLE)/(TOTAL_NUMBER_NODES_XI(2)-1)
318  XI3=0
319  XI3delta=(1.0)/(TOTAL_NUMBER_NODES_XI(3)-1)
320  zero=0
321  DO k=1, TOTAL_NUMBER_NODES_XI(3)
322    !Apex nodes
323    j=1
324    i=1
325    node_idx=node_idx+1
326    CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
327    IF(NodeDomain==ComputationalNodeNumber) THEN
328      FibreFieldAngle=(/zero,zero,zero/) 
329      DO component_idx=1,FieldFibreNumberOfComponents
330        CALL CMISSField_ParameterSetUpdateNode(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
331          & DerivativeUserNumber, &
332          & node_idx,component_idx,FibreFieldAngle(component_idx),Err)
333      ENDDO
334    ENDIF
335    theta=atan(FIBRE_SLOPE_CHANGE*XI3+FIBRE_SLOPE_INTERSECTION)
336    DO j=2, TOTAL_NUMBER_NODES_XI(2) 
337      nu=PI-XI2delta*(j-1)
338      omega=PI/2+cos(2*nu)*atan(SHEET_SLOPE_BASE_ENDO)*(-2*XI3+1)
339      DO i=1, TOTAL_NUMBER_NODES_XI(1)
340        node_idx=node_idx+1
341        CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
342        IF(NodeDomain==ComputationalNodeNumber) THEN
343          FibreFieldAngle=(/theta,zero,omega/)
344          DO component_idx=1,FieldFibreNumberOfComponents
345            CALL CMISSField_ParameterSetUpdateNode(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
346              & 1,DerivativeUserNumber, node_idx,component_idx,FibreFieldAngle(component_idx),Err)
347          ENDDO
348        ENDIF
349      ENDDO
350    ENDDO
351    XI3=XI3+XI3delta
352  ENDDO
353
354  !Create a material field and attach it to the geometric field  
355  CALL CMISSField_Initialise(MaterialField,Err)
356  CALL CMISSField_CreateStart(FieldMaterialUserNumber,Region,MaterialField,Err)
357  CALL CMISSField_TypeSet(MaterialField,CMISS_FIELD_MATERIAL_TYPE,Err)
358  CALL CMISSField_MeshDecompositionSet(MaterialField,Decomposition,Err)        
359  CALL CMISSField_GeometricFieldSet(MaterialField,GeometricField,Err)
360  CALL CMISSField_NumberOfVariablesSet(MaterialField,FieldMaterialNumberOfVariables,Err)
361  CALL CMISSField_NumberOfComponentsSet(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialNumberOfComponents,Err)  
362  CALL CMISSField_ComponentInterpolationSet(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
363  CALL CMISSField_ComponentInterpolationSet(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
364  CALL CMISSField_CreateFinish(MaterialField,Err)
365
366  !Set Mooney-Rivlin constants c10 and c01 to 2.0 and 6.0 respectively.
367  CALL CMISSField_ComponentValuesInitialise(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,2.0_CMISSDP,Err)
368  CALL CMISSField_ComponentValuesInitialise(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,6.0_CMISSDP,Err)
369
370  !Create the equations_set
371  CALL CMISSField_Initialise(EquationsSetField,Err)
372  CALL CMISSEquationsSet_Initialise(EquationsSet,Err)
373  CALL CMISSEquationsSet_CreateStart(EquationSetUserNumber,Region,FibreField,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
374    & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_NO_SUBTYPE,EquationsSetFieldUserNumber,&
375    & EquationsSetField,EquationsSet,Err)
376  ! CMISS_EQUATIONS_SET_ORTHOTROPIC_MATERIAL_COSTA_SUBTYPE
377  CALL CMISSEquationsSet_CreateFinish(EquationsSet,Err)
378
379  !Create the dependent field with 2 variables and 4 components (3 displacement, 1 pressure)
380  CALL CMISSField_Initialise(DependentField,Err)
381  CALL CMISSField_CreateStart(FieldDependentUserNumber,Region,DependentField,Err)
382  CALL CMISSField_TypeSet(DependentField,CMISS_FIELD_GENERAL_TYPE,Err)
383  CALL CMISSField_MeshDecompositionSet(DependentField,Decomposition,Err)
384  CALL CMISSField_GeometricFieldSet(DependentField,GeometricField,Err)
385  CALL CMISSField_DependentTypeSet(DependentField,CMISS_FIELD_DEPENDENT_TYPE,Err)
386  CALL CMISSField_NumberOfVariablesSet(DependentField,FieldDependentNumberOfVariables,Err)
387  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentNumberOfComponents,Err)
388  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,FieldDependentNumberOfComponents,Err)
389  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
390  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
391  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
392  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err)
393  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
394  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
395  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
396  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err)
397  CALL CMISSField_ScalingTypeSet(DependentField,CMISS_FIELD_UNIT_SCALING,Err)
398  CALL CMISSField_CreateFinish(DependentField,Err)
399
400  CALL CMISSEquationsSet_DependentCreateStart(EquationsSet,FieldDependentUserNumber,DependentField,Err)
401  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSet,Err)
402
403  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSet,FieldMaterialUserNumber,MaterialField,Err)  
404  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSet,Err)
405
406  !Create the equations set equations
407  CALL CMISSEquations_Initialise(Equations,Err)
408  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSet,Equations,Err)
409  CALL CMISSEquations_SparsityTypeSet(Equations,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
410  CALL CMISSEquations_OutputTypeSet(Equations,CMISS_EQUATIONS_NO_OUTPUT,Err)
411  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSet,Err)   
412
413  !Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure
414  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
415    & 1,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,Err)
416  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
417    & 2,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,Err)
418  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
419    & 3,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,Err)
420  CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4, &
421    & -14.0_CMISSDP, &
422    & Err)
423
424  !Define the problem
425  CALL CMISSProblem_Initialise(Problem,Err)
426  CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
427  CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_ELASTICITY_CLASS,CMISS_PROBLEM_FINITE_ELASTICITY_TYPE, &
428    & CMISS_PROBLEM_NO_SUBTYPE,Err)
429  CALL CMISSProblem_CreateFinish(Problem,Err)
430
431  !Create the problem control loop
432  CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
433  CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
434
435  !Create the problem solvers
436  CALL CMISSSolver_Initialise(Solver,Err)
437  CALL CMISSSolver_Initialise(LinearSolver,Err)
438  CALL CMISSProblem_SolversCreateStart(Problem,Err)
439  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,1,Solver,Err)
440  CALL CMISSSolver_OutputTypeSet(Solver,CMISS_SOLVER_PROGRESS_OUTPUT,Err)
441  !CALL CMISSSolver_NewtonJacobianCalculationTypeSet(Solver,CMISS_SOLVER_NEWTON_JACOBIAN_FD_CALCULATED,Err)
442  CALL CMISSSolver_NewtonJacobianCalculationTypeSet(Solver,CMISS_SOLVER_NEWTON_JACOBIAN_EQUATIONS_CALCULATED,Err)
443  CALL CMISSSolver_NewtonLinearSolverGet(Solver,LinearSolver,Err)
444  CALL CMISSSolver_LinearTypeSet(LinearSolver,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
445  CALL CMISSProblem_SolversCreateFinish(Problem,Err)
446
447  !Create the problem solver equations
448  CALL CMISSSolver_Initialise(Solver,Err)
449  CALL CMISSSolverEquations_Initialise(SolverEquations,Err)
450  CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)   
451  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,1,Solver,Err)
452  CALL CMISSSolver_SolverEquationsGet(Solver,SolverEquations,Err)
453  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquations,CMISS_SOLVER_SPARSE_MATRICES,Err)
454  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquations,EquationsSet,EquationsSetIndex,Err)
455  CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
456
457  !Prescribe boundary conditions (absolute nodal parameters)
458  CALL CMISSBoundaryConditions_Initialise(BoundaryConditions,Err)
459  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquations,BoundaryConditions,Err)
460
461  !Grab the list of nodes on inner, outer and top surfaces
462  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_TOP_SURFACE,TopSurfaceNodes,TopNormalXi,Err)
463  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_INNER_SURFACE,InnerSurfaceNodes,InnerNormalXi,Err)
464  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_OUTER_SURFACE,OuterSurfaceNodes,OuterNormalXi,Err)
465
466  ! ASSIGN BOUNDARY CONDITIONS
467  !Fix base of the ellipsoid in z direction
468  DO NN=1,SIZE(TopSurfaceNodes,1)
469    NODE=TopSurfaceNodes(NN)
470    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
471    IF(NodeDomain==ComputationalNodeNumber) THEN
472      CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, &
473        & ZCoord, &
474        & Err)
475      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, &
476        & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err)
477    ENDIF
478  ENDDO
479
480  !Apply inner surface pressure
481  !NOTE: Surface pressure goes into pressure_values_set_type of the DELUDELN type
482  DO NN=1,SIZE(InnerSurfaceNodes,1)
483    NODE=InnerSurfaceNodes(NN)
484    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
485    IF(NodeDomain==ComputationalNodeNumber) THEN
486      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, &
487        & ABS(InnerNormalXi), &
488        & CMISS_BOUNDARY_CONDITION_PRESSURE,INNER_PRESSURE,Err)
489    ENDIF
490  ENDDO
491
492  !Apply outer surface pressure
493  DO NN=1,SIZE(OuterSurfaceNodes,1)
494    NODE=OuterSurfaceNodes(NN)
495    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
496    IF(NodeDomain==ComputationalNodeNumber) THEN
497      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, &
498        & ABS(OuterNormalXi), &
499        & CMISS_BOUNDARY_CONDITION_PRESSURE,OUTER_PRESSURE,Err)
500    ENDIF
501  ENDDO
502
503  !Fix more nodes at the base to stop free body motion
504  X_FIXED=.FALSE.
505  Y_FIXED=.FALSE.
506  DO NN=1,SIZE(TopSurfaceNodes,1)
507    NODE=TopSurfaceNodes(NN)
508    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
509    IF(NodeDomain==ComputationalNodeNumber) THEN
510      CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,1, &
511        & XCoord, &
512        & Err)
513      CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,2, &
514        & YCoord, &
515        & Err)
516      IF(ABS(XCoord)<1.0E-6_CMISSDP) THEN
517        CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,1, &
518          & CMISS_BOUNDARY_CONDITION_FIXED,XCoord,Err)
519        WRITE(*,*) "FIXING NODE",NODE,"IN X DIRECTION"
520        X_FIXED=.TRUE.
521      ENDIF
522      IF(ABS(YCoord)<1.0E-6_CMISSDP) THEN
523        CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,2, &
524          & CMISS_BOUNDARY_CONDITION_FIXED,YCoord,Err)
525        WRITE(*,*) "FIXING NODE",NODE,"IN Y DIRECTION"
526        Y_FIXED=.TRUE.
527      ENDIF
528    ENDIF
529  ENDDO
530  CALL MPI_REDUCE(X_FIXED,X_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
531  CALL MPI_REDUCE(Y_FIXED,Y_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
532  IF(ComputationalNodeNumber==0) THEN
533    IF(.NOT.(X_OKAY.AND.Y_OKAY)) THEN
534      WRITE(*,*) "Free body motion could not be prevented!"
535      CALL CMISSFinalise(Err)
536      STOP
537    ENDIF
538  ENDIF
539
540  CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquations,Err)
541
542  !Solve problem
543  CALL CMISSProblem_Solve(Problem,Err)
544
545  !Output solution  
546  CALL CMISSFields_Initialise(Fields,Err)
547  CALL CMISSFields_Create(Region,Fields,Err)
548  CALL CMISSFields_NodesExport(Fields,"QuadraticEllipsoid","FORTRAN",Err)
549  CALL CMISSFields_ElementsExport(Fields,"QuadraticEllipsoid","FORTRAN",Err)
550  CALL CMISSFields_Finalise(Fields,Err)
551
552  CALL CMISSFinalise(Err)
553
554  WRITE(*,'(A)') "Program successfully completed."
555
556  STOP
557
558END PROGRAM QUADRATICELLIPSOIDEEXAMPLE
559