PageRenderTime 130ms CodeModel.GetById 24ms app.highlight 97ms RepoModel.GetById 1ms app.codeStats 0ms

/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/IncompressibleElasticityDrivenMultiCompDarcy/src/IncompressibleElasticityDrivenMultiCompDarcyExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 1658 lines | 896 code | 202 blank | 560 comment | 0 complexity | be21c1a47c3d354ec4cf4a3847da36e9 MD5 | raw file

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

  1!> \file
  2!> \author Christian Michler, Adam Reeve, Andrew Cookson
  3!> \brief This is an example program to solve a coupled Finite Elastiticity Darcy 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):
 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/IncompressibleElasticityDrivenDarcy/src/IncompressibleElasticityDrivenDarcyExample.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/IncompressibleElasticityDrivenDarcy/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/IncompressibleElasticityDrivenDarcy/build-intel'>Linux GNU Build</a>
 47!!
 48!<
 49
 50! !
 51! !  This example considers a coupled Finite Elasticity Multi-Compartment Darcy problem
 52! !
 53
 54!> Main program
 55
 56PROGRAM FINITEELASTICITYMULTICOMPDARCYEXAMPLE
 57
 58  !
 59  !================================================================================================================================
 60  !
 61
 62  !PROGRAM LIBRARIES
 63
 64  USE OPENCMISS
 65  USE FLUID_MECHANICS_IO_ROUTINES
 66  USE MPI
 67
 68#ifdef WIN32
 69  USE IFQWINCMISS
 70#endif
 71
 72  !
 73  !================================================================================================================================
 74  !
 75
 76  !PROGRAM VARIABLES AND TYPES
 77
 78  IMPLICIT NONE
 79
 80  !Test program parameters
 81
 82  REAL(CMISSDP), PARAMETER :: Y_DIM=1.0_CMISSDP
 83  REAL(CMISSDP), PARAMETER :: X_DIM=1.0_CMISSDP
 84  REAL(CMISSDP), PARAMETER :: Z_DIM=3.0_CMISSDP
 85
 86  INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=1
 87  INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=2
 88  INTEGER(CMISSIntg), PARAMETER :: CubicBasisUserNumber=3
 89
 90  INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
 91  INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=2
 92  INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=3
 93  INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=4
 94  INTEGER(CMISSIntg), PARAMETER :: GeometricFieldUserNumber=5
 95  INTEGER(CMISSIntg), PARAMETER :: DependentFieldUserNumberDarcy=6
 96  INTEGER(CMISSIntg) :: MaterialsFieldUserNumberDarcy
 97  INTEGER(CMISSIntg) :: EquationsSetUserNumberDarcy
 98  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=14
 99  INTEGER(CMISSIntg) :: EquationsSetFieldUserNumberDarcy
100  INTEGER(CMISSIntg) :: icompartment,Ncompartments,num_var,componentnum,Nparams
101  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1
102  INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2
103  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1
104  INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1
105  INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1
106  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPorosity=1
107  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPermOverVis=2
108  INTEGER(CMISSIntg) :: SourceFieldDarcyUserNumber
109
110  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1
111  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
112  INTEGER(CMISSIntg) :: IndependentFieldDarcyUserNumber
113  !Program types
114
115  INTEGER(CMISSIntg) :: NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS,NUMBER_GLOBAL_Z_ELEMENTS
116
117!   INTEGER(CMISSIntg) :: MPI_IERROR
118  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,NumberOfDomains,ComputationalNodeNumber
119
120  TYPE(EXPORT_CONTAINER):: CM
121
122  !Program variables
123
124  INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
125  INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
126  INTEGER(CMISSIntg) :: RESTART_VALUE
127
128  INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT
129  INTEGER(CMISSIntg) :: COMPONENT_NUMBER
130  INTEGER(CMISSIntg) :: NODE_NUMBER
131  INTEGER(CMISSIntg) :: ELEMENT_NUMBER
132  INTEGER(CMISSIntg) :: CONDITION
133
134  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY
135  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE
136  INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE
137
138  REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z
139  REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2
140  REAL(CMISSDP) :: GEOMETRY_TOLERANCE
141  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID
142  REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4)
143  REAL(CMISSDP) :: INITIAL_FIELD_SOLID(4)
144  REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
145  REAL(CMISSDP) :: RELATIVE_TOLERANCE
146  REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
147  REAL(CMISSDP) :: LINESEARCH_ALPHA
148  REAL(CMISSDP) :: VALUE
149  REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY
150
151  LOGICAL :: EXPORT_FIELD_IO
152  LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG
153
154  !CMISS variables
155
156  !Regions
157  TYPE(CMISSRegionType) :: Region
158  TYPE(CMISSRegionType) :: WorldRegion
159  !Coordinate systems
160  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem
161  TYPE(CMISSCoordinateSystemType) :: WorldCoordinateSystem
162  !Basis
163  TYPE(CMISSBasisType) :: CubicBasis, QuadraticBasis, LinearBasis, Bases(2)
164  !Meshes
165  TYPE(CMISSMeshType) :: Mesh
166  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
167  !Decompositions
168  TYPE(CMISSDecompositionType) :: Decomposition
169  !Fields
170  TYPE(CMISSFieldsType) :: Fields
171  !Field types
172  TYPE(CMISSFieldType) :: GeometricField
173  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: MaterialsFieldDarcy
174  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: EquationsSetFieldDarcy
175  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: SourceFieldDarcy
176  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: IndependentFieldDarcy
177  !Boundary conditions
178  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy
179  !Equations sets
180  TYPE(CMISSEquationsSetType), ALLOCATABLE, DIMENSION(:) :: EquationsSetDarcy
181  !Equations
182  TYPE(CMISSEquationsType), ALLOCATABLE, DIMENSION(:) :: EquationsDarcy
183  !Problems
184  TYPE(CMISSProblemType) :: Problem
185  !Control loops
186  TYPE(CMISSControlLoopType) :: ControlLoop
187  !Solvers
188  TYPE(CMISSSolverType) :: DynamicSolverDarcy
189  TYPE(CMISSSolverType) :: LinearSolverDarcy
190!   TYPE(CMISSSolverType) :: LinearSolverSolid
191  !Solver equations
192  TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy
193  !Other variables
194  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face1Nodes(:),Face2Nodes(:)
195  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face3Nodes(:),Face4Nodes(:)
196  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face5Nodes(:),Face6Nodes(:)
197  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face7Nodes(:),Face8Nodes(:)
198  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face9Nodes(:),Face10Nodes(:)
199  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face11Nodes(:),Face12Nodes(:)
200  INTEGER(CMISSIntg) :: FaceXi(6)
201  INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
202  REAL(CMISSDP) :: XCoord,YCoord,ZCoord
203  LOGICAL :: X_FIXED,Y_FIXED !,X_OKAY,Y_OKAY
204
205#ifdef WIN32
206  !Quickwin type
207  LOGICAL :: QUICKWIN_STATUS=.FALSE.
208  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
209#endif
210
211  !Generic CMISS variables
212
213  INTEGER(CMISSIntg) :: EquationsSetIndex
214  INTEGER(CMISSIntg) :: Err
215  !Array containing the field variable types that will be used (for ease of incorporating inside a loop)
216  INTEGER(CMISSIntg), ALLOCATABLE, DIMENSION(:) :: VariableTypes
217  REAL(CMISSDP), ALLOCATABLE, DIMENSION(:,:) :: CouplingCoeffs,ConstitutiveParams
218
219  INTEGER(CMISSIntg) :: DIAG_LEVEL_LIST(5)
220!   CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(8) !,TIMING_ROUTINE_LIST(1)
221  CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(1) !,TIMING_ROUTINE_LIST(1)
222
223  !
224  !--------------------------------------------------------------------------------------------------------------------------------
225  !
226
227  !Program variables and types (finite elasticity part)
228
229  !Test program parameters
230
231  INTEGER(CMISSIntg) :: SolidMeshComponenetNumber
232
233  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidUserNumber=51
234  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfVariables=1
235  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfComponents=3
236
237  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidUserNumber=52
238  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfVariables=1
239  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfComponents=3
240
241  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidUserNumber=53
242  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfVariables=1
243  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfComponents=3
244
245  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidUserNumber=54
246  INTEGER(CMISSIntg) :: FieldDependentSolidNumberOfVariables
247  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4
248  INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4  !(u,v,w,m)
249
250  INTEGER(CMISSIntg), PARAMETER :: EquationSetSolidUserNumber=55
251  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldSolidUserNumber=25
252
253  INTEGER(CMISSIntg), PARAMETER :: SolidDisplMeshComponentNumber=1
254  INTEGER(CMISSIntg), PARAMETER :: SolidLagrMultMeshComponentNumber=2
255  INTEGER(CMISSIntg), PARAMETER :: SolidGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
256
257  INTEGER(CMISSIntg), PARAMETER :: DarcyVelMeshComponentNumber=SolidLagrMultMeshComponentNumber
258  INTEGER(CMISSIntg), PARAMETER :: DarcyMassIncreaseMeshComponentNumber=SolidLagrMultMeshComponentNumber
259!   INTEGER(CMISSIntg), PARAMETER :: DarcyGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
260
261  INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=32
262  !Program types
263  !Program variables
264
265  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME
266  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME
267  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA
268  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
269
270  !CMISS variables
271
272!   TYPE(CMISSBasisType) :: BasisSolid
273  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsSolid
274  TYPE(CMISSEquationsType) :: EquationsSolid
275  TYPE(CMISSEquationsSetType) :: EquationsSetSolid
276  TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid,DependentFieldSolid,EquationsSetFieldSolid
277  TYPE(CMISSSolverType) :: SolverSolid
278  TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid
279  TYPE(CMISSMeshElementsType) :: MeshElementsSolid
280
281  !End - Program variables and types (finite elasticity part)
282
283  INTEGER(CMISSIntg) :: dummy_counter
284
285  !
286  !--------------------------------------------------------------------------------------------------------------------------------
287  !
288
289
290#ifdef WIN32
291  !Initialise QuickWin
292  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
293  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
294  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
295  !Set the window parameters
296  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
297  !If attempt fails set with system estimated values
298  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
299#endif
300
301  !
302  !================================================================================================================================
303  !
304
305  !INITIALISE OPENCMISS
306
307  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
308
309  !CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
310
311  !
312  !================================================================================================================================
313  !
314
315  !PROBLEM CONTROL PANEL
316  NUMBER_GLOBAL_X_ELEMENTS=1
317  NUMBER_GLOBAL_Y_ELEMENTS=1
318  NUMBER_GLOBAL_Z_ELEMENTS=3
319
320  IF(NUMBER_GLOBAL_Z_ELEMENTS==0)THEN
321    NUMBER_OF_DIMENSIONS=2
322  ELSE
323    NUMBER_OF_DIMENSIONS=3
324  ENDIF
325  !PROBLEM CONTROL PANEL
326
327!   BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION
328  BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION
329  !Set geometric tolerance
330  GEOMETRY_TOLERANCE = 1.0E-12_CMISSDP
331  !Set initial values
332  INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
333  INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
334  INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
335  INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
336  !Set material parameters
337  POROSITY_PARAM_DARCY=0.1_CMISSDP
338  PERM_OVER_VIS_PARAM_DARCY=1.0_CMISSDP
339  !Set output parameter
340  !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
341  DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_MATRIX_OUTPUT
342  LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_NO_OUTPUT
343  !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
344  EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
345
346  !Set time parameter
347  DYNAMIC_SOLVER_DARCY_START_TIME=1.0E-3_CMISSDP
348  DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP
349  DYNAMIC_SOLVER_DARCY_STOP_TIME=2_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
350  DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
351  !Set result output parameter
352  DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
353  !Set solver parameters
354  LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
355  RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
356  ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
357  DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
358  MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
359  RESTART_VALUE=30_CMISSIntg !default: 30
360  LINESEARCH_ALPHA=1.0_CMISSDP
361
362!   !Import cmHeart mesh information
363!   CALL FLUID_MECHANICS_IO_READ_CMHEART(CM,Err)
364!   BASIS_NUMBER_GEOMETRY=CM%ID_M
365!   BASIS_NUMBER_VELOCITY=CM%ID_V
366!   BASIS_NUMBER_PRESSURE=CM%ID_P
367!   NUMBER_OF_DIMENSIONS=CM%D
368!   BASIS_TYPE=CM%IT_T
369!   BASIS_XI_INTERPOLATION_GEOMETRY=CM%IT_M
370!   BASIS_XI_INTERPOLATION_VELOCITY=CM%IT_V
371!   BASIS_XI_INTERPOLATION_PRESSURE=CM%IT_P
372!   BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION
373!   NUMBER_OF_NODES_GEOMETRY=CM%N_M
374!   NUMBER_OF_NODES_VELOCITY=CM%N_V
375!   NUMBER_OF_NODES_PRESSURE=CM%N_P
376!   TOTAL_NUMBER_OF_NODES=CM%N_T
377!   TOTAL_NUMBER_OF_ELEMENTS=CM%E_T
378!   NUMBER_OF_ELEMENT_NODES_GEOMETRY=CM%EN_M
379!   NUMBER_OF_ELEMENT_NODES_VELOCITY=CM%EN_V
380!   NUMBER_OF_ELEMENT_NODES_PRESSURE=CM%EN_P
381!   !Set domain dimensions
382!   DOMAIN_X1 = -5.0_CMISSDP
383!   DOMAIN_X2 =  5.0_CMISSDP
384!   DOMAIN_Y1 = -5.0_CMISSDP
385!   DOMAIN_Y2 =  5.0_CMISSDP
386!   DOMAIN_Z1 = -5.0_CMISSDP
387!   DOMAIN_Z2 =  5.0_CMISSDP
388!   !Set geometric tolerance
389!   GEOMETRY_TOLERANCE = 1.0E-12_CMISSDP
390!   !Set initial values
391!   INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
392!   INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
393!   INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
394!   INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
395!   INITIAL_FIELD_MAT_PROPERTIES(1)=0.0_CMISSDP
396!   INITIAL_FIELD_MAT_PROPERTIES(2)=0.0_CMISSDP
397!   INITIAL_FIELD_MAT_PROPERTIES(3)=0.0_CMISSDP
398! !   INITIAL_FIELD_SOLID(1)=1.0_CMISSDP
399! !   INITIAL_FIELD_SOLID(2)=1.0_CMISSDP
400! !   INITIAL_FIELD_SOLID(3)=1.0_CMISSDP
401! !   INITIAL_FIELD_SOLID(4)=1.0_CMISSDP
402!   !Set material parameters
403!   POROSITY_PARAM_DARCY=0.1_CMISSDP
404!   PERM_OVER_VIS_PARAM_DARCY=1.0e-1_CMISSDP
405!   POROSITY_PARAM_MAT_PROPERTIES=POROSITY_PARAM_DARCY
406!   PERM_OVER_VIS_PARAM_MAT_PROPERTIES=PERM_OVER_VIS_PARAM_DARCY
407!   !Set number of Gauss points (Mind that also material field may be interpolated)
408!   BASIS_XI_GAUSS_GEOMETRY=3 !4
409!   BASIS_XI_GAUSS_VELOCITY=3 !4
410!   BASIS_XI_GAUSS_PRESSURE=3 !4
411!   !Set output parameter
412!   !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
413!   LINEAR_SOLVER_MAT_PROPERTIES_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
414!   DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
415!   LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_SOLVER_OUTPUT
416!   !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
417!   EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
418!   EQUATIONS_MAT_PROPERTIES_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
419! 
420!   !Set time parameter
421!   DYNAMIC_SOLVER_DARCY_START_TIME=0.0_CMISSDP
422! !   DYNAMIC_SOLVER_DARCY_STOP_TIME=0.03_CMISSDP
423!   DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-2_CMISSDP
424!   DYNAMIC_SOLVER_DARCY_STOP_TIME=1_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
425!   DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
426!   !Set result output parameter
427!   DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
428!   !Set solver parameters
429!   LINEAR_SOLVER_MAT_PROPERTIES_DIRECT_FLAG=.TRUE.
430!   LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
431!   RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
432!   ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
433!   DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
434!   MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
435!   RESTART_VALUE=30_CMISSIntg !default: 30
436!   LINESEARCH_ALPHA=1.0_CMISSDP
437
438  icompartment =1_CMISSIntg
439  Ncompartments=2_CMISSIntg
440  !
441  !================================================================================================================================
442  !
443
444! !   !Set diagnostics
445! ! 
446! !   DIAG_LEVEL_LIST(1)=1
447! !   DIAG_LEVEL_LIST(2)=2
448! !   DIAG_LEVEL_LIST(3)=3
449! !   DIAG_LEVEL_LIST(4)=4
450! !   DIAG_LEVEL_LIST(5)=5
451! ! 
452! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_FINITE_ELEMENT_CALCULATE"
453! ! !   DIAG_ROUTINE_LIST(2)="DARCY_EQUATION_PRE_SOLVE_STORE_REFERENCE_DATA"
454! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_GET_SOLID_DISPLACEMENT"
455! ! !   DIAG_ROUTINE_LIST(2)="DARCY_EQUATION_PRE_SOLVE_ALE_UPDATE_MESH"
456! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_GET_SOLID_DISPLACEMENT"
457! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_UPDATE_BOUNDARY_CONDITIONS"
458! !   DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
459! ! !   DIAG_ROUTINE_LIST(2)="FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR"
460! ! !   DIAG_ROUTINE_LIST(3)="EVALUATE_CHAPELLE_PIOLA_TENSOR_ADDITION"
461! ! !   DIAG_ROUTINE_LIST(5)="DARCY_EQUATION_PRE_SOLVE_MAT_PROPERTIES"
462! ! !   DIAG_ROUTINE_LIST(6)="DATA_FITTING_FINITE_ELEMENT_CALCULATE"
463! ! !   DIAG_ROUTINE_LIST(7)="FINITE_ELASTICITY_FINITE_ELEMENT_JACOBIAN_EVALUATE"
464! ! !   DIAG_ROUTINE_LIST(8)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
465! ! !   DIAG_ROUTINE_LIST(1)="PROBLEM_SOLVER_EQUATIONS_SOLVE"
466! ! !   DIAG_ROUTINE_LIST(1)="SOLVER_NEWTON_SOLVE"
467! ! !   DIAG_ROUTINE_LIST(2)="SOLVER_NEWTON_LINESEARCH_SOLVE"
468! ! !   DIAG_ROUTINE_LIST(1)="SOLVER_SOLUTION_UPDATE"
469! ! !   DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
470
471
472!Set diagnostics
473
474  DIAG_LEVEL_LIST(1)=1
475  DIAG_LEVEL_LIST(2)=2
476  DIAG_LEVEL_LIST(3)=3
477  DIAG_LEVEL_LIST(4)=4
478  DIAG_LEVEL_LIST(5)=5
479
480  !DIAG_ROUTINE_LIST(1)="WRITE_IP_INFO"
481!   DIAG_ROUTINE_LIST(2)="FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR"
482  DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
483
484  !CMISS_ALL_DIAG_TYPE/CMISS_IN_DIAG_TYPE/CMISS_FROM_DIAG_TYPE
485  CALL CMISSDiagnosticsSetOn(CMISS_IN_DIAG_TYPE,DIAG_LEVEL_LIST,"Diagnostics",DIAG_ROUTINE_LIST,Err)
486
487  !CMISS_ALL_TIMING_TYPE/CMISS_IN_TIMING_TYPE/CMISS_FROM_TIMING_TYPE
488  !TIMING_ROUTINE_LIST(1)="PROBLEM_FINITE_ELEMENT_CALCULATE"
489  !CALL TIMING_SET_ON(IN_TIMING_TYPE,.TRUE.,"",TIMING_ROUTINE_LIST,ERR,ERROR,*999)
490
491!   !CMISS_ALL_DIAG_TYPE/CMISS_IN_DIAG_TYPE/CMISS_FROM_DIAG_TYPE
492!    CALL CMISSDiagnosticsSetOn(CMISS_IN_DIAG_TYPE,DIAG_LEVEL_LIST,"Diagnostics",DIAG_ROUTINE_LIST,Err)
493
494  !CMISS_ALL_TIMING_TYPE/CMISS_IN_TIMING_TYPE/CMISS_FROM_TIMING_TYPE
495  !TIMING_ROUTINE_LIST(1)="PROBLEM_FINITE_ELEMENT_CALCULATE"
496  !CALL TIMING_SET_ON(IN_TIMING_TYPE,.TRUE.,"",TIMING_ROUTINE_LIST,ERR,ERROR,*999)
497  ALLOCATE (EquationsSetDarcy(Ncompartments))
498  ALLOCATE (EquationsSetFieldDarcy(Ncompartments))
499  ALLOCATE (MaterialsFieldDarcy(Ncompartments))
500  ALLOCATE (EquationsDarcy(Ncompartments))
501  ALLOCATE (SourceFieldDarcy(Ncompartments))
502  ALLOCATE (IndependentFieldDarcy(Ncompartments))
503  !
504  !================================================================================================================================
505  !
506
507  !Get the number of computational nodes and this computational node number
508  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
509  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
510
511  NumberOfDomains = NumberOfComputationalNodes
512  write(*,*) "NumberOfDomains = ",NumberOfDomains
513
514
515  !
516  !================================================================================================================================
517  !
518
519  !COORDINATE SYSTEM
520
521  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
522  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
523  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NUMBER_OF_DIMENSIONS,Err)
524  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
525
526  !
527  !================================================================================================================================
528  !
529
530  !REGION
531  !For a volume-coupled problem, solid and fluid are based in the same region
532
533  CALL CMISSRegion_Initialise(Region,Err)
534  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
535  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
536  CALL CMISSRegion_CreateFinish(Region,Err)
537
538  !
539  !================================================================================================================================
540  !
541
542  !BASES
543  !Define basis functions
544  CALL CMISSBasis_Initialise(LinearBasis,Err)
545  CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
546  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
547    & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
548  !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err)
549  CALL CMISSBasis_CreateFinish(LinearBasis,Err)
550
551  CALL CMISSBasis_Initialise(QuadraticBasis,Err)
552  CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
553  CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
554    & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
555  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
556    & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
557  !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err)
558  CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
559
560  CALL CMISSBasis_Initialise(CubicBasis,Err)
561  CALL CMISSBasis_CreateStart(CubicBasisUserNumber,CubicBasis,Err)
562  CALL CMISSBasis_InterpolationXiSet(CubicBasis,(/CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION, &
563    & CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION/),Err)
564  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(CubicBasis, &
565    & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
566  !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(CubicBasis,.true.,Err) !Enable 3D interpolation on faces
567  CALL CMISSBasis_CreateFinish(CubicBasis,Err)
568
569  !LinearBasis/QuadraticBasis/CubicBasis
570  Bases(1)=QuadraticBasis
571  Bases(2)=LinearBasis
572!   Bases(1)=CubicBasis
573!   Bases(2)=QuadraticBasis
574
575  !Start the creation of a generated mesh in the region
576  CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
577  CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
578  CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_REGULAR_MESH_TYPE,Err)
579  CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,Bases,Err)
580  IF(NUMBER_GLOBAL_Z_ELEMENTS==0) THEN
581    CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/X_DIM,Y_DIM/),Err)
582    CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS/),Err)
583  ELSE
584    CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/X_DIM,Y_DIM,Z_DIM/),Err)
585    CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS, &
586      & NUMBER_GLOBAL_Z_ELEMENTS/),Err)
587  ENDIF
588  CALL CMISSMesh_Initialise(Mesh,Err)
589  CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
590
591  !
592  !================================================================================================================================
593  !
594
595  !GEOMETRIC FIELD
596
597  !Create a decomposition:
598  !All mesh components (associated with G.Projection / Darcy / solid) share the same decomposition
599  CALL CMISSDecomposition_Initialise(Decomposition,Err)
600  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
601  !Set the decomposition to be a general decomposition with the specified number of domains
602  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
603  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
604  CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.true.,Err)
605  !Finish the decomposition
606  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
607
608  CALL CMISSField_Initialise(GeometricField,Err)
609  CALL CMISSField_CreateStart(GeometricFieldUserNumber,Region,GeometricField,Err)
610  CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
611  CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
612  CALL CMISSField_NumberOfVariablesSet(GeometricField,FieldGeometryNumberOfVariables,Err)
613  CALL CMISSField_NumberOfComponentsSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,Err)  
614  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
615  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
616  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
617  CALL CMISSField_CreateFinish(GeometricField,Err)
618  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricField,Err)
619
620  !--------------------------------------------------------------------------------------------------------------------------------
621  ! Solid
622
623  !Create a decomposition
624
625  !Create a field to put the geometry (defualt is geometry)
626
627  SolidMeshComponenetNumber = SolidGeometryMeshComponentNumber
628
629  CALL CMISSField_Initialise(GeometricFieldSolid,Err)
630  CALL CMISSField_CreateStart(FieldGeometrySolidUserNumber,Region,GeometricFieldSolid,Err)
631  CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err)
632  CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err)
633  CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometrySolidNumberOfVariables,Err)
634  CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometrySolidNumberOfComponents,Err)
635  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidMeshComponenetNumber,Err)
636  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidMeshComponenetNumber,Err)
637  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidMeshComponenetNumber,Err)
638  CALL CMISSField_CreateFinish(GeometricFieldSolid,Err)
639  !Set the mesh component to be used by the field components.
640  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err)
641!   !Set the scaling to use
642!   CALL CMISSField_ScalingTypeSet(GeometricFieldSolid,CMISS_FIELD_NO_SCALING,Err)
643
644  !Create a fibre field and attach it to the geometric field
645  CALL CMISSField_Initialise(FibreFieldSolid,Err)
646  CALL CMISSField_CreateStart(FieldFibreSolidUserNumber,Region,FibreFieldSolid,Err)
647  CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err)
648  CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err)
649  CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err)
650  CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreSolidNumberOfVariables,Err)
651  CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreSolidNumberOfComponents,Err)
652  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
653  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
654  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
655  CALL CMISSField_CreateFinish(FibreFieldSolid,Err)
656
657  ! end Solid
658  !--------------------------------------------------------------------------------------------------------------------------------
659
660  !
661  !================================================================================================================================
662  !
663
664  !EQUATIONS SETS
665  DO icompartment = 1,Ncompartments
666    EquationsSetFieldUserNumberDarcy = 100_CMISSIntg+icompartment
667    EquationsSetUserNumberDarcy = 200_CMISSIntg+icompartment
668  !Create the equations set for ALE Darcy
669    CALL CMISSField_Initialise(EquationsSetFieldDarcy(icompartment),Err)
670    CALL CMISSEquationsSet_Initialise(EquationsSetDarcy(icompartment),Err)
671    CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricField, &
672      & CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
673      & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
674      & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy(icompartment),EquationsSetDarcy(icompartment),Err)
675    !Finish creating the equations set
676    CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy(icompartment),Err)
677    !Set the values for the equations set field to be the current compartment number (1 - N), and the total number of compartments (N)
678    CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
679      & CMISS_FIELD_VALUES_SET_TYPE,1,icompartment,Err)
680    CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
681      & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
682  ENDDO
683
684  !--------------------------------------------------------------------------------------------------------------------------------
685  ! Solid
686
687  !Create the equations_set
688  CALL CMISSField_Initialise(EquationsSetFieldSolid,Err)
689  CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err)
690  CALL CMISSEquationsSet_CreateStart(EquationSetSolidUserNumber,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
691    & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
692    & EquationsSetFieldSolidUserNumber,EquationsSetFieldSolid,EquationsSetSolid,Err)
693  CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err)
694  !Set the values for the equations set field to be the current compartment number (O for the finite elasticity equations_set), and the total number of compartments (N)
695  !Need to store number of compartments, as finite elasticity uses this to calculate the total mass increase for the constiutive law
696  CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
697     & CMISS_FIELD_VALUES_SET_TYPE,1,0_CMISSIntg,Err)
698  CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
699     & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
700
701
702  ! end Solid
703  !--------------------------------------------------------------------------------------------------------------------------------
704  !--------------------------------------------------------------------------------------------------------------------------------
705  ! Solid Materials Field
706
707  !Create a material field and attach it to the geometric field
708  CALL CMISSField_Initialise(MaterialFieldSolid,Err)
709  !
710  CALL CMISSField_CreateStart(FieldMaterialSolidUserNumber,Region,MaterialFieldSolid,Err)
711  !
712  CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err)
713  CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err)
714  CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err)
715  CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialSolidNumberOfVariables,Err)
716  CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialSolidNumberOfComponents,Err)
717  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
718  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
719  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
720  !
721  CALL CMISSField_CreateFinish(MaterialFieldSolid,Err)
722
723  !Set material parameters
724  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
725    & 2.0_CMISSDP,Err)
726!   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,2.0e3_CMISSDP,Err)
727  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2, &
728    & 6.0_CMISSDP,Err)
729!   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,33.0_CMISSDP,Err)
730  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3, &
731    & 10.0_CMISSDP,Err)
732
733  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialSolidUserNumber,MaterialFieldSolid,Err)
734  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err)
735
736  ! end Solid
737  !--------------------------------------------------------------------------------------------------------------------------------
738
739
740  !--------------------------------------------------------------------------------------------------------------------------------
741  ! Solid
742
743  !Create a dependent field with two variables and four components
744  CALL CMISSField_Initialise(DependentFieldSolid,Err)
745  !
746  CALL CMISSField_CreateStart(FieldDependentSolidUserNumber,Region,DependentFieldSolid,Err)
747  !
748  CALL CMISSField_TypeSet(DependentFieldSolid,CMISS_FIELD_GENERAL_TYPE,Err)
749  CALL CMISSField_MeshDecompositionSet(DependentFieldSolid,Decomposition,Err)
750  CALL CMISSField_GeometricFieldSet(DependentFieldSolid,GeometricFieldSolid,Err)
751  CALL CMISSField_DependentTypeSet(DependentFieldSolid,CMISS_FIELD_DEPENDENT_TYPE,Err)
752  !Create 2N+2 number of variables - 2 for solid, 2N for N Darcy compartments
753  FieldDependentSolidNumberOfVariables=2*Ncompartments+2
754  CALL CMISSField_NumberOfVariablesSet(DependentFieldSolid,FieldDependentSolidNumberOfVariables,Err)
755  !create two variables for each compartment
756  ALLOCATE(VariableTypes(2*Ncompartments+2))
757  DO num_var=1,Ncompartments+1
758     VariableTypes(2*num_var-1)=CMISS_FIELD_U_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
759     VariableTypes(2*num_var)=CMISS_FIELD_DELUDELN_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
760  ENDDO
761  CALL CMISSField_VariableTypesSet(DependentFieldSolid,VariableTypes,Err) 
762!   CALL CMISSField_VariableTypesSet(DependentFieldSolid,(/CMISS_FIELD_U_VARIABLE_TYPE, &
763!     & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err)
764    CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
765       & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
766    CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
767       & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
768  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
769  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
770    & FieldDependentSolidNumberOfComponents,Err)
771!   DO icompartment=3,2*Ncompartments+2
772!     CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
773!   ENDDO
774!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
775!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
776!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
777
778  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidDisplMeshComponentNumber,Err)
779  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidDisplMeshComponentNumber,Err)
780  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidDisplMeshComponentNumber,Err)
781  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4, &
782    & CMISS_FIELD_NODE_BASED_INTERPOLATION, &
783    & Err)
784!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
785!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
786  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidLagrMultMeshComponentNumber,Err)
787!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
788!     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
789!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
790!     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
791!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
792!     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
793  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
794    & SolidDisplMeshComponentNumber, &
795    & Err)
796  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
797    & SolidDisplMeshComponentNumber, &
798    & Err)
799  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
800    & SolidDisplMeshComponentNumber, &
801    & Err)
802  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
803    & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
804!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
805!     & CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
806!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
807  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
808    & SolidLagrMultMeshComponentNumber, &
809    & Err)
810  !loop over the number of compartments
811  DO icompartment=3,2*Ncompartments+2
812!     CALL CMISSField_DimensionSet(DependentFieldSolid,VariableTypes(icompartment), &
813!        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
814    CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
815    DO componentnum=1,FieldDependentFluidNumberOfComponents-1
816    !set dimension type
817!     CALL CMISSField_DimensionSet(DependentField,VariableTypes(icompartment), &
818!        & CMISS_FIELD_SCALAR_DIMENSION_TYPE,Err)
819      CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, &
820       & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
821      CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, & 
822         & DarcyVelMeshComponentNumber,Err)
823    ENDDO
824      CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment), &
825       & FieldDependentFluidNumberOfComponents, &
826       & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
827!     CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
828!       & FieldDependentFluidNumberOfComponents,MESH_COMPONENT_NUMBER_PRESSURE,Err)
829    CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
830      & FieldDependentFluidNumberOfComponents,DarcyMassIncreaseMeshComponentNumber,Err)
831    
832  ENDDO
833
834!   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
835!   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
836!   !For this equation type, MESH_COMPONENT_NUMBER_PRESSURE is actually the mass increase component as the pressure is taken from the solid equations
837!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
838!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
839!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
840!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
841!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
842!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
843!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
844!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
845
846  !
847  CALL CMISSField_CreateFinish(DependentFieldSolid,Err)
848  !
849  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentSolidUserNumber,DependentFieldSolid,Err)
850  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err)
851
852!   !Initialise dependent field (solid displacement and pressure)
853!   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS !+1
854!     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
855!       & COMPONENT_NUMBER,INITIAL_FIELD_SOLID(COMPONENT_NUMBER),Err)
856!   ENDDO
857
858  ! end Solid
859  !--------------------------------------------------------------------------------------------------------------------------------
860
861  DO icompartment = 1,Ncompartments
862    CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy(icompartment),FieldDependentSolidUserNumber,&
863      & DependentFieldSolid,Err)
864    CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy(icompartment),Err)
865  ENDDO
866
867
868  !Create the equations set dependent field variables for ALE Darcy
869!   CALL CMISSField_Initialise(DependentFieldDarcy,Err)
870!   CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,DependentFieldUserNumberDarcy, & ! ??? UserNumber ???
871!   CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,FieldDependentSolidUserNumber, & ! ??? UserNumber ???
872!     & DependentFieldSolid,Err)
873! !   !Set the mesh component to be used by the field components.
874! !   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
875! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
876! !       & MESH_COMPONENT_NUMBER_VELOCITY,Err)
877! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
878! !       & MESH_COMPONENT_NUMBER_VELOCITY,Err)
879! !   ENDDO
880! !   COMPONENT_NUMBER=NUMBER_OF_DIMENSIONS+1
881! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
882! !       & MESH_COMPONENT_NUMBER_PRESSURE,Err)
883! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
884! !       & MESH_COMPONENT_NUMBER_PRESSURE,Err)
885!   !Finish the equations set dependent field variables
886!   CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy,Err)
887
888  !Initialise dependent field (velocity components,pressure,mass increase)
889  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS+1
890    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
891      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
892    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
893      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
894  ENDDO
895  !
896  !================================================================================================================================
897  !
898  ALLOCATE(CouplingCoeffs(Ncompartments,Ncompartments))
899  IF(Ncompartments==2)THEN
900    CouplingCoeffs(1,1)=0.0E-01_CMISSDP
901!     CouplingCoeffs(1,2)=-1.0E-04_CMISSDP
902!     CouplingCoeffs(2,1)=-1.0E-04_CMISSDP
903    CouplingCoeffs(1,2)=0.0E-01_CMISSDP
904    CouplingCoeffs(2,1)=0.0E-01_CMISSDP
905    CouplingCoeffs(2,2)=0.0E-01_CMISSDP
906  ELSE IF(Ncompartments==3)THEN
907    CouplingCoeffs(1,1)=1.0E-02_CMISSDP
908    CouplingCoeffs(1,2)=1.0E-02_CMISSDP
909    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
910    CouplingCoeffs(2,1)=1.0E-02_CMISSDP
911    CouplingCoeffs(2,2)=2.0E-02_CMISSDP
912    CouplingCoeffs(2,3)=1.0E-02_CMISSDP
913    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
914    CouplingCoeffs(3,2)=1.0E-02_CMISSDP
915    CouplingCoeffs(3,3)=1.0E-02_CMISSDP
916  ELSE IF(Ncompartments==4)THEN
917    CouplingCoeffs(1,1)=0.0E-02_CMISSDP
918    CouplingCoeffs(1,2)=0.0E-02_CMISSDP
919    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
920    CouplingCoeffs(1,4)=0.0E-02_CMISSDP
921    CouplingCoeffs(2,1)=0.0E-02_CMISSDP
922    CouplingCoeffs(2,2)=0.0E-02_CMISSDP
923    CouplingCoeffs(2,3)=0.0E-02_CMISSDP
924    CouplingCoeffs(2,4)=0.0E-02_CMISSDP
925    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
926    CouplingCoeffs(3,2)=0.0E-02_CMISSDP
927    CouplingCoeffs(3,3)=0.0E-02_CMISSDP
928    CouplingCoeffs(3,4)=0.0E-02_CMISSDP
929    CouplingCoeffs(4,1)=0.0E-02_CMISSDP
930    CouplingCoeffs(4,2)=0.0E-02_CMISSDP
931    CouplingCoeffs(4,3)=0.0E-02_CMISSDP
932    CouplingCoeffs(4,4)=0.0E-02_CMISSDP
933  ELSE IF(Ncompartments==5)THEN
934    CouplingCoeffs(1,1)=0.0E-02_CMISSDP
935    CouplingCoeffs(1,2)=0.0E-02_CMISSDP
936    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
937    CouplingCoeffs(1,4)=0.0E-02_CMISSDP
938    CouplingCoeffs(1,5)=0.0E-02_CMISSDP
939    CouplingCoeffs(2,1)=0.0E-02_CMISSDP
940    CouplingCoeffs(2,2)=0.0E-02_CMISSDP
941    CouplingCoeffs(2,3)=0.0E-02_CMISSDP
942    CouplingCoeffs(2,4)=0.0E-02_CMISSDP
943    CouplingCoeffs(2,5)=0.0E-02_CMISSDP
944    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
945    CouplingCoeffs(3,2)=0.0E-02_CMISSDP
946    CouplingCoeffs(3,3)=0.0E-02_CMISSDP
947    CouplingCoeffs(3,4)=0.0E-02_CMISSDP
948    CouplingCoeffs(3,5)=0.0E-02_CMISSDP
949    CouplingCoeffs(4,1)=0.0E-02_CMISSDP
950    CouplingCoeffs(4,2)=0.0E-02_CMISSDP
951    CouplingCoeffs(4,3)=0.0E-02_CMISSDP
952    CouplingCoeffs(4,4)=0.0E-02_CMISSDP
953    CouplingCoeffs(4,5)=0.0E-02_CMISSDP
954    CouplingCoeffs(5,1)=0.0E-02_CMISSDP
955    CouplingCoeffs(5,2)=0.0E-02_CMISSDP
956    CouplingCoeffs(5,3)=0.0E-02_CMISSDP
957    CouplingCoeffs(5,4)=0.0E-02_CMISSDP
958    CouplingCoeffs(5,5)=0.0E-02_CMISSDP
959  ELSE
960    write(*,*) "Can't initialise coupling coefficients array."
961  ENDIF
962  !Define the material parameters for each compartments' constitutive law (for determining pressure)
963  Nparams=3
964  ALLOCATE(ConstitutiveParams(Ncompartments,Nparams))
965  IF(Ncompartments==2)THEN
966    ConstitutiveParams(1,1)=10.0E-01_CMISSDP
967    ConstitutiveParams(1,2)=10.0E-01_CMISSDP
968    ConstitutiveParams(1,3)=10.0E-01_CMISSDP
969    ConstitutiveParams(2,1)=10.0E-01_CMISSDP
970    ConstitutiveParams(2,2)=10.0E-01_CMISSDP
971    ConstitutiveParams(2,3)=10.0E-01_CMISSDP
972  ELSE IF(Ncompartments==3)THEN
973    ConstitutiveParams(1,1)=1.0E-02_CMISSDP
974    ConstitutiveParams(1,2)=1.0E-02_CMISSDP
975    ConstitutiveParams(1,3)=0.0E-02_CMISSDP
976    ConstitutiveParams(2,1)=1.0E-02_CMISSDP
977    ConstitutiveParams(2,2)=2.0E-02_CMISSDP
978    ConstitutiveParams(2,3)=1.0E-02_CMISSDP
979    ConstitutiveParams(3,1)=0.0E-02_CMISSDP
980    ConstitutiveParams(3,2)=1.0E-02_CMISSDP
981    ConstitutiveParams(3,3)=1.0E-02_CMISSDP
982  ELSE IF(Ncompartments==4)THEN
983    ConstitutiveParams(1,1)=0.0E-02_CMISSDP
984    ConstitutiveParams(1,2)=0.0E-02_CMISSDP
985    ConstitutiveParams(1,3)=0.0E-02_CMISSDP
986    ConstitutiveP

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