PageRenderTime 110ms CodeModel.GetById 17ms app.highlight 83ms RepoModel.GetById 3ms app.codeStats 0ms

/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenDarcy/src/QuadraticEllipsoidDrivenDarcyExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 1124 lines | 649 code | 181 blank | 294 comment | 0 complexity | 065c8a8617ed85fdc58477d6f0c89bff MD5 | raw file

Large files files are truncated, but you can click here to view the full file

  1!> \file
  2!> \author Christian Michler
  3!> \brief This is an example program to solve a coupled Finite Elastiticity Darcy equation on a cylindrical geometry 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): Christian Michler, 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 MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenDarcy/src/QuadraticEllipsoidDrivenDarcyExample.f90
 43!! Example program to solve coupled FiniteElasticityDarcy equations using OpenCMISS calls.
 44!! \par Latest Builds:
 45!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenDarcy/build-intel'>Linux Intel Build</a>
 46!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenDarcy/build-intel'>Linux GNU Build</a>
 47!!
 48!<
 49
 50! !
 51! !  This example considers a coupled Finite Elasticity Darcy problem on a cylindrical geometry
 52! !
 53
 54!> Main program
 55
 56PROGRAM QUADRATICELLIPSOIDDRIVENDARCYEXAMPLE
 57
 58  !
 59  !================================================================================================================================
 60  !
 61
 62  !PROGRAM LIBRARIES
 63
 64  USE OPENCMISS
 65  USE MPI
 66
 67#ifdef WIN32
 68  USE IFQWINCMISS
 69#endif
 70
 71  !
 72  !================================================================================================================================
 73  !
 74
 75  !PROGRAM VARIABLES AND TYPES
 76
 77  IMPLICIT NONE
 78
 79  !Test program parameters
 80
 81  REAL(CMISSDP), PARAMETER :: PI=3.14159_CMISSDP
 82  REAL(CMISSDP), PARAMETER :: LONG_AXIS=2.0_CMISSDP
 83  REAL(CMISSDP), PARAMETER :: SHORT_AXIS=1.0_CMISSDP
 84  REAL(CMISSDP), PARAMETER :: WALL_THICKNESS=0.5_CMISSDP
 85  REAL(CMISSDP), PARAMETER :: CUTOFF_ANGLE=1.5708_CMISSDP
 86  REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_INTERSECTION=1.73205_CMISSDP !Slope of fibres in base endocardium = 60 degrees
 87  REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_CHANGE=-3.4641_CMISSDP !Slope change of fibres from 60 to -60 degrees in transmural direction 
 88  REAL(CMISSDP), PARAMETER :: SHEET_SLOPE_BASE_ENDO=1.0_CMISSDP !Slope of sheet at base endocardium 
 89
 90  REAL(CMISSDP), PARAMETER :: INNER_PRESSURE=2.0_CMISSDP  !Positive is compressive
 91  REAL(CMISSDP), PARAMETER :: OUTER_PRESSURE=0.0_CMISSDP !Positive is compressive
 92  REAL(CMISSDP), PARAMETER :: C1= 2.0_CMISSDP
 93  REAL(CMISSDP), PARAMETER :: C2= 6.0_CMISSDP
 94  REAL(CMISSDP), PARAMETER :: C3=10.0_CMISSDP
 95
 96  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalXElements=4  ! X ==NUMBER_GLOBAL_CIRCUMFERENTIAL_ELEMENTS
 97  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalYElements=4  ! Y ==NUMBER_GLOBAL_LONGITUDINAL_ELEMENTS
 98  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalZElements=2  ! Z ==NUMBER_GLOBAL_TRANSMURAL_ELEMENTS
 99  INTEGER(CMISSIntg) :: NumberOfDomains
100
101  INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
102  INTEGER(CMISSIntg), PARAMETER :: NumberOfSpatialCoordinates=3
103  INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=1
104  INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=1
105  INTEGER(CMISSIntg), PARAMETER :: QuadraticCollapsedBasisUserNumber=2
106  INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=3
107  INTEGER(CMISSIntg), PARAMETER :: LinearCollapsedBasisUserNumber=4
108  INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=1
109  INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=2
110  INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=1
111  INTEGER(CMISSIntg), PARAMETER :: DerivativeUserNumber=1
112
113  INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshDimensions=3
114  INTEGER(CMISSIntg), PARAMETER :: NumberOfXiCoordinates=3
115  INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshComponents=2
116  INTEGER(CMISSIntg), PARAMETER :: QuadraticMeshComponentNumber=1
117  INTEGER(CMISSIntg), PARAMETER :: LinearMeshComponentNumber=2
118  INTEGER(CMISSIntg), PARAMETER :: TotalNumberOfElements=1
119
120  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumberSolid=1
121  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumberDarcy=11
122  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1
123  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
124
125  INTEGER(CMISSIntg), PARAMETER :: FieldFibreUserNumber=2
126  INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfVariables=1
127  INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfComponents=3
128
129  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialUserNumber=3
130  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfVariables=1
131  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfComponents=3 !2
132
133  INTEGER(CMISSIntg), PARAMETER :: FieldDependentUserNumber=4
134  INTEGER(CMISSIntg), PARAMETER :: FieldDependentNumberOfVariables=4 !2
135  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4
136  INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4
137
138  INTEGER(CMISSIntg), PARAMETER :: IndependentFieldDarcyUserNumber=15
139
140  INTEGER(CMISSIntg), PARAMETER :: EquationSetUserNumberSolid=1
141  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumberSolid=13
142  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=1
143
144
145  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcy=8
146  INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberDarcy=12
147  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumberDarcy=22
148  INTEGER(CMISSIntg), PARAMETER :: SourceFieldDarcyUserNumber=42
149
150  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1
151  INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2
152  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1
153  INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1
154  INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1
155
156  !Program types
157
158
159  !Program variables
160
161  INTEGER(CMISSIntg) :: MPI_IERROR
162  INTEGER(CMISSIntg) :: EquationsSetIndex  
163  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber
164  REAL(CMISSDP) :: FibreFieldAngle(3) 
165  REAL(CMISSDP) :: nu,theta,omega,XI3,XI3delta,XI2delta, zero
166  INTEGER(CMISSIntg) ::i,j,k,component_idx,node_idx,TOTAL_NUMBER_NODES_XI(3)
167  !For grabbing surfaces
168  INTEGER(CMISSIntg) :: InnerNormalXi,OuterNormalXi,TopNormalXi
169  INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodes(:)
170  INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodes(:)
171  INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodes(:)
172  INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodesDarcyVel(:)
173  INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodesDarcyVel(:)
174  INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodesDarcyVel(:)
175  INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
176  REAL(CMISSDP) :: XCoord,YCoord,ZCoord
177  LOGICAL :: X_FIXED,Y_FIXED,X_OKAY,Y_OKAY
178
179  INTEGER(CMISSIntg) :: GeometricFieldDarcyMeshComponentNumber, DarcyVelMeshComponentNumber, DarcyMassIncreaseMeshComponentNumber
180  INTEGER(CMISSIntg) :: MeshComponentNumber_dummy
181
182  INTEGER(CMISSIntg) :: NUMBER_OF_DOMAINS
183
184  INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
185
186  INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
187  INTEGER(CMISSIntg) :: RESTART_VALUE
188
189  INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT
190  INTEGER(CMISSIntg) :: COMPONENT_NUMBER, NODE_NUMBER
191  INTEGER(CMISSIntg) :: CONDITION
192
193  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY
194  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE
195  INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE
196
197  REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z
198  REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2
199  REAL(CMISSDP) :: GEOMETRY_TOLERANCE
200  INTEGER(CMISSIntg) :: EDGE_COUNT
201  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID
202  REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4)
203  REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
204  REAL(CMISSDP) :: RELATIVE_TOLERANCE
205  REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
206  REAL(CMISSDP) :: LINESEARCH_ALPHA
207  REAL(CMISSDP) :: VALUE
208  REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY
209
210  LOGICAL :: EXPORT_FIELD_IO
211  LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG
212
213  !CMISS variables
214
215  TYPE(CMISSBasisType) :: QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis
216  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditions
217  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem, WorldCoordinateSystem
218  TYPE(CMISSMeshType) :: Mesh
219  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
220  TYPE(CMISSDecompositionType) :: Decomposition
221  TYPE(CMISSEquationsType) :: Equations
222  TYPE(CMISSEquationsSetType) :: EquationsSetSolid
223  TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid
224  TYPE(CMISSFieldType) :: DependentField,EquationsSetFieldSolid
225
226  TYPE(CMISSFieldType) :: IndependentFieldDarcy
227
228  TYPE(CMISSFieldType) :: GeometricFieldDarcy
229
230  TYPE(CMISSFieldType) :: MaterialsFieldDarcy
231  TYPE(CMISSFieldType) :: EquationsSetFieldDarcy
232  TYPE(CMISSFieldType) :: SourceFieldDarcy
233  !Boundary conditions
234  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy
235  !Equations sets
236  TYPE(CMISSEquationsSetType) :: EquationsSetDarcy
237  !Equations
238  TYPE(CMISSEquationsType) :: EquationsDarcy
239
240  TYPE(CMISSFieldsType) :: Fields
241  TYPE(CMISSProblemType) :: Problem
242  TYPE(CMISSRegionType) :: Region,WorldRegion
243  TYPE(CMISSSolverType) :: SolverSolid,LinearSolverSolid
244  TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid
245  TYPE(CMISSControlLoopType) :: ControlLoop
246
247  !Solvers
248  TYPE(CMISSSolverType) :: DynamicSolverDarcy
249  TYPE(CMISSSolverType) :: LinearSolverDarcy
250  !Solver equations
251  TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy
252
253  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME
254  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME
255  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA
256  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
257
258  INTEGER(CMISSIntg),allocatable :: ElementUserNodes(:)
259  INTEGER(CMISSIntg) :: NUMBER_USER_ELEMENT_NODES, ELEMENT_NUMBER
260
261#ifdef WIN32
262  !Quickwin type
263  LOGICAL :: QUICKWIN_STATUS=.FALSE.
264  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
265#endif
266
267  !Generic CMISS variables
268  INTEGER(CMISSIntg) :: Err
269
270#ifdef WIN32
271  !Initialise QuickWin
272  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
273  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
274  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
275  !Set the window parameters
276  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
277  !If attempt fails set with system estimated values
278  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
279#endif
280
281  !Set initial values
282  INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
283  INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
284  INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
285  INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
286  !Set material parameters
287  POROSITY_PARAM_DARCY=0.1_CMISSDP
288  PERM_OVER_VIS_PARAM_DARCY=1.0_CMISSDP
289  !Set output parameter
290  !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
291  DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
292  LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_SOLVER_OUTPUT
293  !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
294  EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
295
296  !Set time parameter
297  DYNAMIC_SOLVER_DARCY_START_TIME=1.0E-3_CMISSDP
298  DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP
299  DYNAMIC_SOLVER_DARCY_STOP_TIME=2_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
300  DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
301  !Set result output parameter
302  DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
303  !Set solver parameters
304  LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
305  RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
306  ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
307  DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
308  MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
309  RESTART_VALUE=30_CMISSIntg !default: 30
310  LINESEARCH_ALPHA=1.0_CMISSDP
311
312
313  !LinearMeshComponentNumber/QuadraticMeshComponentNumber
314  DarcyVelMeshComponentNumber = LinearMeshComponentNumber
315  DarcyMassIncreaseMeshComponentNumber = LinearMeshComponentNumber
316!   GeometricFieldDarcyMeshComponentNumber = DarcyVelMeshComponentNumber
317  GeometricFieldDarcyMeshComponentNumber = QuadraticMeshComponentNumber
318
319
320  !
321  !================================================================================================================================
322  !
323
324  !Intialise cmiss
325  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
326
327  CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
328
329  WRITE(*,'(A)') "Program starting."
330
331  !Set all diganostic levels on for testing
332  CALL CMISSDiagnosticsSetOn(CMISS_FROM_DIAG_TYPE,(/1,2,3,4,5/),"Diagnostics",(/"PROBLEM_FINITE_ELEMENT_CALCULATE"/),Err)
333
334  !Get the number of computational nodes and this computational node number
335  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
336  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
337
338  !Broadcast the number of elements in the X,Y and Z directions and the number of partitions to the other computational nodes
339!   CALL MPI_BCAST(NumberGlobalXElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
340!   CALL MPI_BCAST(NumberGlobalYElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
341!   CALL MPI_BCAST(NumberGlobalZElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
342!   CALL MPI_BCAST(NumberOfDomains,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
343  NumberOfDomains=NumberOfComputationalNodes
344
345  !Create a CS - default is 3D rectangular cartesian CS with 0,0,0 as origin
346  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
347  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
348  CALL CMISSCoordinateSystem_TypeSet(CoordinateSystem,CMISS_COORDINATE_RECTANGULAR_CARTESIAN_TYPE,Err)
349  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NumberOfSpatialCoordinates,Err)
350  CALL CMISSCoordinateSystem_OriginSet(CoordinateSystem,(/0.0_CMISSDP,0.0_CMISSDP,0.0_CMISSDP/),Err)
351  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
352
353  !
354  !================================================================================================================================
355  !
356
357  !Create a region and assign the CS to the region
358  CALL CMISSRegion_Initialise(Region,Err)
359  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
360  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
361  CALL CMISSRegion_CreateFinish(Region,Err)
362
363  !
364  !================================================================================================================================
365  !
366
367  !Define basis functions - tri-linear Lagrange and tri-Quadratic Lagrange, each with collapsed variant
368    !Quadratic Basis
369  CALL CMISSBasis_Initialise(QuadraticBasis,Err)
370  CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
371  CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
372    & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
373  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
374    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
375!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
376  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err) !Have to do this
377  CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
378  
379    !Collapsed Quadratic Basis
380  CALL CMISSBasis_Initialise(QuadraticCollapsedBasis,Err)
381  CALL CMISSBasis_CreateStart(QuadraticCollapsedBasisUserNumber,QuadraticCollapsedBasis,Err)
382  CALL CMISSBasis_TypeSet(QuadraticCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
383  CALL CMISSBasis_NumberOfXiSet(QuadraticCollapsedBasis,NumberOfXiCoordinates,Err)
384  CALL CMISSBasis_InterpolationXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
385       & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
386  CALL CMISSBasis_CollapsedXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED, &
387       & CMISS_BASIS_COLLAPSED_AT_XI0,CMISS_BASIS_NOT_COLLAPSED/),Err)
388  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticCollapsedBasis, &
389       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)  
390!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
391  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticCollapsedBasis,.true.,Err) !Have to do this
392  CALL CMISSBasis_CreateFinish(QuadraticCollapsedBasis,Err)
393
394    !Linear Basis
395  CALL CMISSBasis_Initialise(LinearBasis,Err)
396  CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
397  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
398    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
399!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
400  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
401  CALL CMISSBasis_CreateFinish(LinearBasis,Err)
402
403    !Collapsed Linear Basis
404  CALL CMISSBasis_Initialise(LinearCollapsedBasis,Err)
405  CALL CMISSBasis_CreateStart(LinearCollapsedBasisUserNumber,LinearCollapsedBasis,Err)
406  CALL CMISSBasis_TypeSet(LinearCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
407  CALL CMISSBasis_NumberOfXiSet(LinearCollapsedBasis,NumberOfXiCoordinates,Err)
408  CALL CMISSBasis_InterpolationXiSet(LinearCollapsedBasis,(/CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION, &
409       & CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION,CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION/),Err)
410  CALL CMISSBasis_CollapsedXiSet(LinearCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED,CMISS_BASIS_COLLAPSED_AT_XI0, &
411    & CMISS_BASIS_NOT_COLLAPSED/),Err)
412  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearCollapsedBasis, &
413       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
414!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
415  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearCollapsedBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
416  CALL CMISSBasis_CreateFinish(LinearCollapsedBasis,Err)
417
418  !
419  !================================================================================================================================
420  !
421
422  !Start the creation of a generated ellipsoid mesh
423  CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
424  CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
425  !Set up an ellipsoid mesh
426  CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_MESH_TYPE,Err)
427  !Set the quadratic and linear bases
428  CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,[QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis],Err)
429  !Define the mesh on the region
430  CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/LONG_AXIS,SHORT_AXIS,WALL_THICKNESS,CUTOFF_ANGLE/),Err)
431  CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NumberGlobalXElements,NumberGlobalYElements, &
432    & NumberGlobalZElements/),Err)
433  
434  !Finish the creation of a generated mesh in the region
435  CALL CMISSMesh_Initialise(Mesh,Err)
436  CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
437
438  !
439  !================================================================================================================================
440  !
441
442  !Create a decomposition
443  CALL CMISSDecomposition_Initialise(Decomposition,Err)
444  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
445  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
446  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
447  CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.TRUE.,Err)
448  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
449
450  !
451  !================================================================================================================================
452  !
453
454  ! --- GeometricFieldSolid ---
455  !Create a field to put the geometry (default is geometry)
456  CALL CMISSField_Initialise(GeometricFieldSolid,Err)
457  CALL CMISSField_CreateStart(FieldGeometryUserNumberSolid,Region,GeometricFieldSolid,Err)
458  CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err)
459  CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err)
460  CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometryNumberOfVariables,Err)
461  CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)  
462  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
463  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
464  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
465  CALL CMISSField_CreateFinish(GeometricFieldSolid,Err)
466
467  !Update the geometric field parameters
468  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err)
469
470  !
471  !================================================================================================================================
472  !
473
474  ! --- GeometricFieldDarcy ---
475  !Create a field to put the geometry (default is geometry)
476  CALL CMISSField_Initialise(GeometricFieldDarcy,Err)
477  CALL CMISSField_CreateStart(FieldGeometryUserNumberDarcy,Region,GeometricFieldDarcy,Err)
478  CALL CMISSField_MeshDecompositionSet(GeometricFieldDarcy,Decomposition,Err)
479  CALL CMISSField_TypeSet(GeometricFieldDarcy,CMISS_FIELD_GEOMETRIC_TYPE,Err)
480  CALL CMISSField_NumberOfVariablesSet(GeometricFieldDarcy,FieldGeometryNumberOfVariables,Err)
481  CALL CMISSField_NumberOfComponentsSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)  
482  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,1, &
483    & GeometricFieldDarcyMeshComponentNumber,Err)
484  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,2, &
485    & GeometricFieldDarcyMeshComponentNumber,Err)
486  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,3, &
487    & GeometricFieldDarcyMeshComponentNumber,Err)
488  CALL CMISSField_CreateFinish(GeometricFieldDarcy,Err)
489
490  !Update the geometric field parameters
491  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldDarcy,Err)
492
493  !
494  !================================================================================================================================
495  !
496
497  !Create a fibre field and attach it to the geometric field
498  CALL CMISSField_Initialise(FibreFieldSolid,Err)
499  CALL CMISSField_CreateStart(FieldFibreUserNumber,Region,FibreFieldSolid,Err)
500  CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err)
501  CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err)
502  CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err)
503  CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreNumberOfVariables,Err)
504  CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreNumberOfComponents,Err)  
505  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
506  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
507  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
508  CALL CMISSField_CreateFinish(FibreFieldSolid,Err)
509
510  !Set Fibre directions (this block is parallel-untested)
511  node_idx=0  
512  !This is valid only for quadratic basis functions
513  TOTAL_NUMBER_NODES_XI(1)=NumberGlobalXElements*2
514  TOTAL_NUMBER_NODES_XI(2)=NumberGlobalYElements*2+1
515  TOTAL_NUMBER_NODES_XI(3)=NumberGlobalZElements*2+1
516
517  XI2delta=(PI-CUTOFF_ANGLE)/(TOTAL_NUMBER_NODES_XI(2)-1)
518  XI3=0
519  XI3delta=(1.0)/(TOTAL_NUMBER_NODES_XI(3)-1)
520  zero=0
521  DO k=1, TOTAL_NUMBER_NODES_XI(3)
522    !Apex nodes
523    j=1
524    i=1
525    node_idx=node_idx+1
526    CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
527    IF(NodeDomain==ComputationalNodeNumber) THEN
528      FibreFieldAngle=(/zero,zero,zero/) 
529      DO component_idx=1,FieldFibreNumberOfComponents
530        CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
531          & DerivativeUserNumber,node_idx,component_idx,FibreFieldAngle(component_idx),Err)
532      ENDDO
533    ENDIF
534    theta=atan(FIBRE_SLOPE_CHANGE*XI3+FIBRE_SLOPE_INTERSECTION)
535    DO j=2, TOTAL_NUMBER_NODES_XI(2) 
536      nu=PI-XI2delta*(j-1)
537      omega=PI/2+cos(2*nu)*atan(SHEET_SLOPE_BASE_ENDO)*(-2*XI3+1)
538      DO i=1, TOTAL_NUMBER_NODES_XI(1)
539        node_idx=node_idx+1
540        CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
541        IF(NodeDomain==ComputationalNodeNumber) THEN
542          FibreFieldAngle=(/theta,zero,omega/)
543          DO component_idx=1,FieldFibreNumberOfComponents
544            CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
545              & DerivativeUserNumber, node_idx,component_idx,FibreFieldAngle(component_idx),Err)
546          ENDDO
547        ENDIF
548      ENDDO
549    ENDDO
550    XI3=XI3+XI3delta
551  ENDDO
552
553  !Create a material field and attach it to the geometric field
554  CALL CMISSField_Initialise(MaterialFieldSolid,Err)
555  CALL CMISSField_CreateStart(FieldMaterialUserNumber,Region,MaterialFieldSolid,Err)
556  CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err)
557  CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err)        
558  CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err)
559  CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialNumberOfVariables,Err)
560  CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialNumberOfComponents,Err)  
561  CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
562  CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
563  CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
564  CALL CMISSField_CreateFinish(MaterialFieldSolid,Err)
565
566  !Set Mooney-Rivlin constants 
567  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,C1,Err)
568  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,C2,Err)
569  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,C3,Err)
570
571  !
572  !================================================================================================================================
573  !
574
575  !EQUATIONS SETS
576
577  !Create the equations set for ALE Darcy
578  CALL CMISSField_Initialise(EquationsSetFieldDarcy,Err)
579  CALL CMISSEquationsSet_Initialise(EquationsSetDarcy,Err)
580  CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricFieldDarcy, &
581    & CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
582    & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,&
583    & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy,EquationsSetDarcy,Err)
584  CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy,Err)
585
586  !Create the equations set for the solid
587  CALL CMISSField_Initialise(EquationsSetFieldSolid,Err)
588  CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err)
589  CALL CMISSEquationsSet_CreateStart(EquationSetUserNumberSolid,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
590    & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE, &
591    & EquationsSetFieldUserNumberSolid,EquationsSetFieldSolid,EquationsSetSolid,Err)
592  CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err)
593
594  !
595  !================================================================================================================================
596  !
597
598  !DEPENDENT FIELDS
599
600  !Create a dependent field with four variables (U, DelUDelN = solid, V, DelVDelN = Darcy) and four components
601  !Solid: The U, DelUDelN variables have 4 components (3 displacement, 1 pressure)
602  CALL CMISSField_Initialise(DependentField,Err)
603  CALL CMISSField_CreateStart(FieldDependentUserNumber,Region,DependentField,Err)
604  CALL CMISSField_TypeSet(DependentField,CMISS_FIELD_GENERAL_TYPE,Err)
605  CALL CMISSField_MeshDecompositionSet(DependentField,Decomposition,Err)
606  CALL CMISSField_GeometricFieldSet(DependentField,GeometricFieldSolid,Err)
607  CALL CMISSField_DependentTypeSet(DependentField,CMISS_FIELD_DEPENDENT_TYPE,Err)
608  CALL CMISSField_NumberOfVariablesSet(DependentField,FieldDependentNumberOfVariables,Err)
609
610  CALL CMISSField_VariableTypesSet(DependentField,(/CMISS_FIELD_U_VARIABLE_TYPE, &
611    & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err)
612  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
613  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
614  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
615  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
616
617  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
618  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
619  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
620  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err)
621  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
622  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
623  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
624  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err)
625
626  !Darcy: The V, DelVDelN variables have 4 components (3 velocities, 1 mass increase)
627  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber,Err)
628  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber,Err)
629  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber,Err)
630  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,4,DarcyMassIncreaseMeshComponentNumber,Err)
631  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber,Err)
632  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber,Err)
633  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber,Err)
634  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4, &
635    & DarcyMassIncreaseMeshComponentNumber,Err)
636
637  CALL CMISSField_ScalingTypeSet(DependentField,CMISS_FIELD_UNIT_SCALING,Err)
638
639  CALL CMISSField_CreateFinish(DependentField,Err)
640
641  !
642  !================================================================================================================================
643  !
644
645  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentUserNumber,DependentField,Err)
646  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err)
647
648  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialUserNumber,MaterialFieldSolid,Err)  
649  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err)
650
651  !
652  !================================================================================================================================
653  !
654
655  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,FieldDependentUserNumber,DependentField,Err)
656  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy,Err)
657
658  DO COMPONENT_NUMBER=1,FieldDependentFluidNumberOfComponents
659    CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
660      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
661  ENDDO
662
663  !
664  !================================================================================================================================
665  !
666
667  !INDEPENDENT FIELD Darcy for storing BC flags
668
669  CALL CMISSField_Initialise(IndependentFieldDarcy,Err)
670  CALL CMISSEquationsSet_IndependentCreateStart(EquationsSetDarcy,IndependentFieldDarcyUserNumber, &
671    & IndependentFieldDarcy,Err)
672
673  CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber,Err)
674  CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber,Err)
675  CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber,Err)
676!   CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,4,DarcyMassIncreaseMeshComponentNumber,Err)
677
678  CALL CMISSEquationsSet_IndependentCreateFinish(EquationsSetDarcy,Err)
679
680  CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
681    & 0.0_CMISSDP,Err)
682  CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2, &
683    & 0.0_CMISSDP,Err)
684  CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3, &
685    & 0.0_CMISSDP,Err)
686
687  !
688  !================================================================================================================================
689  !
690
691  !Create the equations set materials field variables for ALE Darcy
692  CALL CMISSField_Initialise(MaterialsFieldDarcy,Err)
693  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy,MaterialsFieldUserNumberDarcy, &
694    & MaterialsFieldDarcy,Err)
695  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy,Err)
696  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
697    & 1,POROSITY_PARAM_DARCY,Err)
698  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
699    & 2,PERM_OVER_VIS_PARAM_DARCY,Err)
700
701  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
702    & 3,0.0_CMISSDP,Err)
703  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
704    & 4,0.0_CMISSDP,Err)
705  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
706    & 5,PERM_OVER_VIS_PARAM_DARCY,Err)
707  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
708    & 6,0.0_CMISSDP,Err)
709  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
710    & 7,PERM_OVER_VIS_PARAM_DARCY,Err)
711
712  !
713  !================================================================================================================================
714  !
715
716  !Source field
717  CALL CMISSField_Initialise(SourceFieldDarcy,Err)
718  CALL CMISSEquationsSet_SourceCreateStart(EquationsSetDarcy,SourceFieldDarcyUserNumber,SourceFieldDarcy,Err)
719  CALL CMISSEquationsSet_SourceCreateFinish(EquationsSetDarcy,Err)
720
721!   ELEMENT_NUMBER = 5
722!   COMPONENT_NUMBER = 4
723!   VALUE = 4.2_CMISSDP
724! !   CALL CMISSField_ParameterSetUpdateElement(RegionUserNumber,SourceFieldDarcyUserNumber,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
725! !     & ELEMENT_NUMBER,COMPONENT_NUMBER,VALUE,Err)
726! 
727!   NUMBER_USER_ELEMENT_NODES = 27  !hardcoding is bad - but how to access the number of element nodes ?
728!                                   !- there is no CMISS library call available for this
729!                                   !- traversing the structure of 'CMISSMeshElementsType' does not work either,
730!                                   !  since certain members are private
731! 
732!   allocate( ElementUserNodes(NUMBER_USER_ELEMENT_NODES) )
733! 
734!   CALL CMISSMeshElements_NodesGet(RegionUserNumber,MeshUserNumber,QuadraticMeshComponentNumber,ELEMENT_NUMBER, &
735!     & ElementUserNodes,Err)
736! 
737!   DO NN=1,NUMBER_USER_ELEMENT_NODES
738!     NODE_NUMBER = ElementUserNodes(NN)
739!     CALL CMISSField_ParameterSetUpdateNode(SourceFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
740!       & 1,NODE_NUMBER,COMPONENT_NUMBER,VALUE,Err)
741!   ENDDO
742
743  !
744  !================================================================================================================================
745  !
746
747  !EQUATIONS SET EQUATIONS
748
749  !Darcy
750  CALL CMISSEquations_Initialise(EquationsDarcy,Err)
751  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy,EquationsDarcy,Err)
752  CALL CMISSEquations_SparsityTypeSet(EquationsDarcy,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
753  CALL CMISSEquations_OutputTypeSet(EquationsDarcy,EQUATIONS_DARCY_OUTPUT,Err)
754  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetDarcy,Err)
755
756  !Solid
757  CALL CMISSEquations_Initialise(Equations,Err)
758  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetSolid,Equations,Err)
759  CALL CMISSEquations_SparsityTypeSet(Equations,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
760  CALL CMISSEquations_OutputTypeSet(Equations,CMISS_EQUATIONS_NO_OUTPUT,Err)
761  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetSolid,Err)
762
763  !
764  !================================================================================================================================
765  !
766
767  !Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure
768  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
769    & CMISS_FIELD_VALUES_SET_TYPE, &
770    & 1,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,Err)
771  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
772    & CMISS_FIELD_VALUES_SET_TYPE, &
773    & 2,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,Err)
774  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
775    & CMISS_FIELD_VALUES_SET_TYPE, &
776    & 3,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,Err)
777!   CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4,-14.0_CMISSDP,Err)
778  CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4,0.0_CMISSDP, &
779    & Err)
780
781  !
782  !================================================================================================================================
783  !
784
785  !Define the problem
786  CALL CMISSProblem_Initialise(Problem,Err)
787  CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
788  CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_MULTI_PHYSICS_CLASS,CMISS_PROBLEM_FINITE_ELASTICITY_DARCY_TYPE, &
789    & CMISS_PROBLEM_QUASISTATIC_ELASTICITY_TRANSIENT_DARCY_SUBTYPE,Err)
790  CALL CMISSProblem_CreateFinish(Problem,Err)
791
792  !
793  !================================================================================================================================
794  !
795
796  !Create the problem control loop
797  CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
798  CALL CMISSControlLoop_Initialise(ControlLoop,Err)
799  CALL CMISSProblem_ControlLoopGet(Problem,CMISS_CONTROL_LOOP_NODE,ControlLoop,Err)
800!   CALL CMISSControlLoop_MaximumIterationsSet(ControlLoop,1,Err)  ! this one sets the increment loop counter
801  CALL CMISSControlLoop_TimesSet(ControlLoop,DYNAMIC_SOLVER_DARCY_START_TIME,DYNAMIC_SOLVER_DARCY_STOP_TIME, &
802    & DYNAMIC_SOLVER_DARCY_TIME_INCREMENT,Err)
803  CALL CMISSControlLoop_TimeOutputSet(ControlLoop,DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY,Err)
804!   CALL CMISSControlLoop_OutputTypeSet(ControlLoop,CMISS_CONTROL_LOOP_PROGRESS_OUTPUT,Err)
805  CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
806
807  !
808  !================================================================================================================================
809  !
810
811  !Create the problem solvers
812  CALL CMISSSolver_Initialise(SolverSolid,Err)
813  CALL CMISSSolver_Initialise(LinearSolverSolid,Err)
814  CALL CMISSSolver_Initialise(DynamicSolverDarcy,Err)
815  CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
816
817  CALL CMISSProblem_SolversCreateStart(Problem,Err)
818
819  ! Solid
820  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
821    & SolverSolidIndex,SolverSolid,Err)
822  CALL CMISSSolver_OutputTypeSet(SolverSolid,CMISS_SOLVER_PROGRESS_OUTPUT,Err)
823!   CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_FD_CALCULATED,Err)
824  CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_EQUATIONS_CALCULATED,Err)
825
826!   CALL CMISSSolverNonLinearTypeSet(SolverSolid,CMISS_SOLVER_NONLINEAR_NEWTON,Err)
827!   CALL CMISSSolver_LibraryTypeSet(SolverSolid,CMISS_SOLVER_PETSC_LIBRARY,Err)
828
829  CALL CMISSSolver_NewtonAbsoluteToleranceSet(SolverSolid,ABSOLUTE_TOLERANCE,Err)
830  CALL CMISSSolver_NewtonRelativeToleranceSet(SolverSolid,RELATIVE_TOLERANCE,Err)
831  CALL CMISSSolver_NewtonMaximumIterationsSet(SolverSolid,MAXIMUM_ITERATIONS,Err)
832
833  CALL CMISSSolver_NewtonLinearSolverGet(SolverSolid,LinearSolverSolid,Err)
834  CALL CMISSSolver_LinearTypeSet(LinearSolverSolid,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
835
836  !Darcy
837  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
838    & SolverDarcyIndex,DynamicSolverDarcy,Err)
839  CALL CMISSSolver_OutputTypeSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE,Err)
840  CALL CMISSSolver_DynamicThetaSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_THETA,Err)
841!   CALL CMISSSolverDynamicDynamicSet(DynamicSolverDarcy,.TRUE.,Err)
842  CALL CMISSSolver_DynamicLinearSolverGet(DynamicSolverDarcy,LinearSolverDarcy,Err)
843  IF(LINEAR_SOLVER_DARCY_DIRECT_FLAG) THEN
844    CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
845    CALL CMISSSolver_LibraryTypeSet(LinearSolverDarcy,CMISS_SOLVER_MUMPS_LIBRARY,Err)
846  ELSE
847    CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
848    CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverDarcy,MAXIMUM_ITERATIONS,Err)
849    CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverDarcy,DIVERGENCE_TOLERANCE,Err)
850    CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverDarcy,RELATIVE_TOLERANCE,Err)
851    CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverDarcy,ABSOLUTE_TOLERANCE,Err)
852    CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverDarcy,RESTART_VALUE,Err)
853  ENDIF
854
855CALL CMISSProblem_SolversCreateFinish(Problem,Err)
856
857
858  !
859  !================================================================================================================================
860  !
861
862  !Create the problem solver equations
863  CALL CMISSSolver_Initialise(SolverSolid,Err)
864  CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
865
866  CALL CMISSSolverEquations_Initialise(SolverEquationsSolid,Err)
867  CALL CMISSSolverEquations_Initialise(SolverEquationsDarcy,Err)
868
869  CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
870  !
871  !Get the finite elasticity solver equations
872  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
873    & SolverSolidIndex,SolverSolid,Err)
874  CALL CMISSSolver_SolverEquationsGet(SolverSolid,SolverEquationsSolid,Err)
875  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsSolid,CMISS_SOLVER_SPARSE_MATRICES,Err)
876  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsSolid,EquationsSetSolid,EquationsSetIndex,Err)
877  !
878  !Get the Darcy solver equations
879  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
880    & SolverDarcyIndex,LinearSolverDarcy,Err)
881  CALL CMISSSolver_SolverEquationsGet(LinearSolverDarcy,SolverEquationsDarcy,Err)
882  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsDarcy,CMISS_SOLVER_SPARSE_MATRICES,Err)
883  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsDarcy,EquationsSetDarcy,EquationsSetIndex,Err)
884  !
885  CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
886
887  !
888  !================================================================================================================================
889  !
890
891  !Prescribe boundary conditions (absolute nodal parameters)
892  CALL CMISSBoundaryConditions_Initialise(BoundaryConditions,Err)
893  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsSolid,BoundaryConditions,Err)
894
895  !Grab the list of nodes on inner, outer and top surfaces
896  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_TOP_SURFACE,TopSurfaceNodes,TopNormalXi,Err)
897  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_INNER_SURFACE,InnerSurfaceNodes,InnerNormalXi,Err)
898  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_OUTER_SURFACE,OuterSurfaceNodes,OuterNormalXi,Err)
899
900  write(*,*)'TopSurfaceNodes = ',TopSurfaceNodes
901
902  ! ASSIGN BOUNDARY CONDITIONS
903  !Fix base of the ellipsoid in z direction
904  DO NN=1,SIZE(TopSurfaceNodes,1)
905    NODE=TopSurfaceNodes(NN)
906    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
907    IF(NodeDomain==ComputationalNodeNumber) THEN
908      CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, &
909        & ZCoord,Err)
910      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, &
911        & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err)
912    ENDIF
913  ENDDO
914
915  !Apply inner surface pressure
916  !NOTE: Surface pressure goes into pressure_values_set_type of the DELUDELN type
917  DO NN=1,SIZE(InnerSurfaceNodes,1)
918    NODE=InnerSurfaceNodes(NN)
919    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
920    IF(NodeDomain==ComputationalNodeNumber) THEN
921      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, &
922        & ABS(InnerNormalXi),CMISS_BOUNDARY_CONDITION_PRESSURE,INNER_PRESSURE,Err)
923    ENDIF
924  ENDDO
925
926  !Apply outer surface pressure
927  DO NN=1,SIZE(OuterSurfaceNodes,1)
928    NODE=OuterSurfaceNodes(NN)
929    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
930    IF(NodeDomain==ComputationalNodeNumber) THEN
931      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, &
932        & ABS(OuterNor

Large files files are truncated, but you can click here to view the full file