PageRenderTime 285ms CodeModel.GetById 78ms app.highlight 190ms RepoModel.GetById 4ms app.codeStats 1ms

/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenMultiCompDarcy/src/QuadraticEllipsoidDrivenMultiCompDarcyExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 1318 lines | 824 code | 154 blank | 340 comment | 0 complexity | 61ea8a05237d93b94245a785e439c9ef 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 QUADRATICELLIPSOIDDRIVENMULTICOMPDARCYEXAMPLE
 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=1  ! 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 :: FieldDependentSolidUserNumber=4
134  INTEGER(CMISSIntg) :: FieldDependentSolidNumberOfVariables
135  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4
136  INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4
137
138  INTEGER(CMISSIntg), PARAMETER :: EquationSetSolidUserNumber=55
139  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldSolidUserNumber=25
140  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=1
141  INTEGER(CMISSIntg) :: EquationsSetFieldUserNumberDarcy
142  INTEGER(CMISSIntg) :: icompartment,Ncompartments,num_var,componentnum,Nparams
143
144  INTEGER(CMISSIntg) :: MaterialsFieldUserNumberDarcy
145  INTEGER(CMISSIntg) :: EquationsSetUserNumberDarcy
146
147
148  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1
149  INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2
150  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1
151  INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1
152  INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1
153  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPorosity=1
154  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPermOverVis=2
155  INTEGER(CMISSIntg), PARAMETER :: SolidDisplMeshComponentNumber=1
156  INTEGER(CMISSIntg), PARAMETER :: SolidLagrMultMeshComponentNumber=2
157  INTEGER(CMISSIntg), PARAMETER :: SolidGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
158  !Program types
159
160
161  !Program variables
162
163  INTEGER(CMISSIntg) :: MPI_IERROR
164  INTEGER(CMISSIntg) :: EquationsSetIndex  
165  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber
166  REAL(CMISSDP) :: FibreFieldAngle(3) 
167  REAL(CMISSDP) :: nu,theta,omega,XI3,XI3delta,XI2delta, zero
168  INTEGER(CMISSIntg) ::i,j,k,component_idx,node_idx,TOTAL_NUMBER_NODES_XI(3)
169  !For grabbing surfaces
170  INTEGER(CMISSIntg) :: InnerNormalXi,OuterNormalXi,TopNormalXi
171  INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodes(:)
172  INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodes(:)
173  INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodes(:)
174  INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodesDarcyVel(:)
175  INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodesDarcyVel(:)
176  INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodesDarcyVel(:)
177  INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
178  REAL(CMISSDP) :: XCoord,YCoord,ZCoord
179  LOGICAL :: X_FIXED,Y_FIXED,X_OKAY,Y_OKAY
180
181  INTEGER(CMISSIntg) :: GeometricFieldDarcyMeshComponentNumber, DarcyVelMeshComponentNumber, DarcyMassIncreaseMeshComponentNumber
182  INTEGER(CMISSIntg) :: MeshComponentNumber_dummy
183
184  INTEGER(CMISSIntg) :: NUMBER_OF_DOMAINS
185
186  INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
187
188  INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
189  INTEGER(CMISSIntg) :: RESTART_VALUE
190
191  INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT
192  INTEGER(CMISSIntg) :: COMPONENT_NUMBER
193  INTEGER(CMISSIntg) :: CONDITION
194
195  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY
196  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE
197  INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE
198
199  REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z
200  REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2
201  REAL(CMISSDP) :: GEOMETRY_TOLERANCE
202  INTEGER(CMISSIntg) :: EDGE_COUNT
203  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID
204  REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4)
205  REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
206  REAL(CMISSDP) :: RELATIVE_TOLERANCE
207  REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
208  REAL(CMISSDP) :: LINESEARCH_ALPHA
209  REAL(CMISSDP) :: VALUE
210  REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY
211
212  LOGICAL :: EXPORT_FIELD_IO
213  LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG
214
215  !CMISS variables
216
217  TYPE(CMISSBasisType) :: QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis
218  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditions
219  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem, WorldCoordinateSystem
220  TYPE(CMISSMeshType) :: Mesh
221  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
222  TYPE(CMISSDecompositionType) :: Decomposition
223  TYPE(CMISSEquationsType) :: Equations
224  TYPE(CMISSEquationsSetType) :: EquationsSetSolid
225  TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid
226  TYPE(CMISSFieldType) :: DependentFieldSolid,EquationsSetFieldSolid
227
228  TYPE(CMISSFieldType) :: GeometricFieldDarcy
229
230  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: MaterialsFieldDarcy
231  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: EquationsSetFieldDarcy
232  !Boundary conditions
233  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy
234  !Equations sets
235  TYPE(CMISSEquationsSetType), ALLOCATABLE, DIMENSION(:) :: EquationsSetDarcy
236  !Equations
237  TYPE(CMISSEquationsType), ALLOCATABLE, DIMENSION(:) :: EquationsDarcy
238
239
240
241  TYPE(CMISSFieldsType) :: Fields
242  TYPE(CMISSProblemType) :: Problem
243  TYPE(CMISSRegionType) :: Region,WorldRegion
244  TYPE(CMISSSolverType) :: SolverSolid,LinearSolverSolid
245  TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid
246  TYPE(CMISSControlLoopType) :: ControlLoop
247
248  !Solvers
249  TYPE(CMISSSolverType) :: DynamicSolverDarcy
250  TYPE(CMISSSolverType) :: LinearSolverDarcy
251  !Solver equations
252  TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy
253
254  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME
255  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME
256  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA
257  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
258
259#ifdef WIN32
260  !Quickwin type
261  LOGICAL :: QUICKWIN_STATUS=.FALSE.
262  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
263#endif
264
265  !Generic CMISS variables
266  INTEGER(CMISSIntg) :: Err
267  INTEGER(CMISSIntg), ALLOCATABLE, DIMENSION(:) :: VariableTypes
268  REAL(CMISSDP), ALLOCATABLE, DIMENSION(:,:) :: CouplingCoeffs,ConstitutiveParams
269#ifdef WIN32
270  !Initialise QuickWin
271  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
272  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
273  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
274  !Set the window parameters
275  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
276  !If attempt fails set with system estimated values
277  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
278#endif
279
280  !Set initial values
281  NUMBER_OF_DIMENSIONS=3_CMISSIntg
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=0.0_CMISSDP
298  DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP
299  DYNAMIC_SOLVER_DARCY_STOP_TIME=5_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  icompartment =1_CMISSIntg
312  Ncompartments=4_CMISSIntg
313
314  !LinearMeshComponentNumber/QuadraticMeshComponentNumber
315  DarcyVelMeshComponentNumber = LinearMeshComponentNumber
316  DarcyMassIncreaseMeshComponentNumber = LinearMeshComponentNumber
317!   GeometricFieldDarcyMeshComponentNumber = DarcyVelMeshComponentNumber
318  GeometricFieldDarcyMeshComponentNumber = QuadraticMeshComponentNumber
319
320  ALLOCATE (EquationsSetDarcy(Ncompartments))
321  ALLOCATE (EquationsSetFieldDarcy(Ncompartments))
322  ALLOCATE (MaterialsFieldDarcy(Ncompartments))
323  ALLOCATE (EquationsDarcy(Ncompartments))
324  !
325  !================================================================================================================================
326  !
327
328  !Intialise cmiss
329  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
330
331  CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
332
333  WRITE(*,'(A)') "Program starting."
334
335  !Set all diganostic levels on for testing
336  CALL CMISSDiagnosticsSetOn(CMISS_FROM_DIAG_TYPE,(/1,2,3,4,5/),"Diagnostics",(/"PROBLEM_FINITE_ELEMENT_CALCULATE"/),Err)
337
338  !Get the number of computational nodes and this computational node number
339  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
340  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
341
342  !Broadcast the number of elements in the X,Y and Z directions and the number of partitions to the other computational nodes
343!   CALL MPI_BCAST(NumberGlobalXElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
344!   CALL MPI_BCAST(NumberGlobalYElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
345!   CALL MPI_BCAST(NumberGlobalZElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
346!   CALL MPI_BCAST(NumberOfDomains,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
347  NumberOfDomains=NumberOfComputationalNodes
348
349  !Create a CS - default is 3D rectangular cartesian CS with 0,0,0 as origin
350  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
351  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
352  CALL CMISSCoordinateSystem_TypeSet(CoordinateSystem,CMISS_COORDINATE_RECTANGULAR_CARTESIAN_TYPE,Err)
353  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NumberOfSpatialCoordinates,Err)
354  CALL CMISSCoordinateSystem_OriginSet(CoordinateSystem,(/0.0_CMISSDP,0.0_CMISSDP,0.0_CMISSDP/),Err)
355  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
356
357  !
358  !================================================================================================================================
359  !
360
361  !Create a region and assign the CS to the region
362  CALL CMISSRegion_Initialise(Region,Err)
363  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
364  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
365  CALL CMISSRegion_CreateFinish(Region,Err)
366
367  !
368  !================================================================================================================================
369  !
370
371  !Define basis functions - tri-linear Lagrange and tri-Quadratic Lagrange, each with collapsed variant
372    !Quadratic Basis
373  CALL CMISSBasis_Initialise(QuadraticBasis,Err)
374  CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
375  CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
376    & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
377  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
378    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
379!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
380  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err) !Have to do this
381  CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
382
383    !Collapsed Quadratic Basis
384  CALL CMISSBasis_Initialise(QuadraticCollapsedBasis,Err)
385  CALL CMISSBasis_CreateStart(QuadraticCollapsedBasisUserNumber,QuadraticCollapsedBasis,Err)
386  CALL CMISSBasis_TypeSet(QuadraticCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
387  CALL CMISSBasis_NumberOfXiSet(QuadraticCollapsedBasis,NumberOfXiCoordinates,Err)
388  CALL CMISSBasis_InterpolationXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
389       & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
390  CALL CMISSBasis_CollapsedXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED, &
391       & CMISS_BASIS_COLLAPSED_AT_XI0,CMISS_BASIS_NOT_COLLAPSED/),Err)
392  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticCollapsedBasis, &
393       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)  
394!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
395  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticCollapsedBasis,.true.,Err) !Have to do this
396  CALL CMISSBasis_CreateFinish(QuadraticCollapsedBasis,Err)
397
398    !Linear Basis
399  CALL CMISSBasis_Initialise(LinearBasis,Err)
400  CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
401  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
402    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
403!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
404  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
405  CALL CMISSBasis_CreateFinish(LinearBasis,Err)
406
407    !Collapsed Linear Basis
408  CALL CMISSBasis_Initialise(LinearCollapsedBasis,Err)
409  CALL CMISSBasis_CreateStart(LinearCollapsedBasisUserNumber,LinearCollapsedBasis,Err)
410  CALL CMISSBasis_TypeSet(LinearCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
411  CALL CMISSBasis_NumberOfXiSet(LinearCollapsedBasis,NumberOfXiCoordinates,Err)
412  CALL CMISSBasis_InterpolationXiSet(LinearCollapsedBasis,(/CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION, &
413       & CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION,CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION/),Err)
414  CALL CMISSBasis_CollapsedXiSet(LinearCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED,CMISS_BASIS_COLLAPSED_AT_XI0, &
415    & CMISS_BASIS_NOT_COLLAPSED/),Err)
416  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearCollapsedBasis, &
417       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
418!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
419  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearCollapsedBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
420  CALL CMISSBasis_CreateFinish(LinearCollapsedBasis,Err)
421
422  !
423  !================================================================================================================================
424  !
425
426  !Start the creation of a generated ellipsoid mesh
427  CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
428  CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
429  !Set up an ellipsoid mesh
430  CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_MESH_TYPE,Err)
431  !Set the quadratic and linear bases
432  CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,[QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis],Err)
433  !Define the mesh on the region
434  CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/LONG_AXIS,SHORT_AXIS,WALL_THICKNESS,CUTOFF_ANGLE/),Err)
435  CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NumberGlobalXElements,NumberGlobalYElements, &
436    & NumberGlobalZElements/),Err)
437  
438  !Finish the creation of a generated mesh in the region
439  CALL CMISSMesh_Initialise(Mesh,Err)
440  CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
441
442  !
443  !================================================================================================================================
444  !
445
446  !Create a decomposition
447  CALL CMISSDecomposition_Initialise(Decomposition,Err)
448  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
449  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
450  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
451  CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.TRUE.,Err)
452  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
453
454  !
455  !================================================================================================================================
456  !
457
458  ! --- GeometricFieldSolid ---
459  !Create a field to put the geometry (default is geometry)
460  CALL CMISSField_Initialise(GeometricFieldSolid,Err)
461  CALL CMISSField_CreateStart(FieldGeometryUserNumberSolid,Region,GeometricFieldSolid,Err)
462  CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err)
463  CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err)
464  CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometryNumberOfVariables,Err)
465  CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)  
466  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
467  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
468  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
469  CALL CMISSField_CreateFinish(GeometricFieldSolid,Err)
470
471  !Update the geometric field parameters
472  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err)
473
474  !
475  !================================================================================================================================
476  !
477
478  ! --- GeometricFieldDarcy ---
479  !Create a field to put the geometry (default is geometry)
480  CALL CMISSField_Initialise(GeometricFieldDarcy,Err)
481  CALL CMISSField_CreateStart(FieldGeometryUserNumberDarcy,Region,GeometricFieldDarcy,Err)
482  CALL CMISSField_MeshDecompositionSet(GeometricFieldDarcy,Decomposition,Err)
483  CALL CMISSField_TypeSet(GeometricFieldDarcy,CMISS_FIELD_GEOMETRIC_TYPE,Err)
484  CALL CMISSField_NumberOfVariablesSet(GeometricFieldDarcy,FieldGeometryNumberOfVariables,Err)
485  CALL CMISSField_NumberOfComponentsSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)  
486  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,1, &
487    & GeometricFieldDarcyMeshComponentNumber,Err)
488  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,2, &
489    & GeometricFieldDarcyMeshComponentNumber,Err)
490  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,3, &
491    & GeometricFieldDarcyMeshComponentNumber,Err)
492  CALL CMISSField_CreateFinish(GeometricFieldDarcy,Err)
493
494  !Update the geometric field parameters
495  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldDarcy,Err)
496
497  !
498  !================================================================================================================================
499  !
500
501  !Create a fibre field and attach it to the geometric field
502  CALL CMISSField_Initialise(FibreFieldSolid,Err)
503  CALL CMISSField_CreateStart(FieldFibreUserNumber,Region,FibreFieldSolid,Err)
504  CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err)
505  CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err)        
506  CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err)
507  CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreNumberOfVariables,Err)
508  CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreNumberOfComponents,Err)  
509  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
510  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
511  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
512  CALL CMISSField_CreateFinish(FibreFieldSolid,Err)
513
514  !Set Fibre directions (this block is parallel-untested)
515  node_idx=0  
516  !This is valid only for quadratic basis functions
517  TOTAL_NUMBER_NODES_XI(1)=NumberGlobalXElements*2
518  TOTAL_NUMBER_NODES_XI(2)=NumberGlobalYElements*2+1
519  TOTAL_NUMBER_NODES_XI(3)=NumberGlobalZElements*2+1
520
521  XI2delta=(PI-CUTOFF_ANGLE)/(TOTAL_NUMBER_NODES_XI(2)-1)
522  XI3=0
523  XI3delta=(1.0)/(TOTAL_NUMBER_NODES_XI(3)-1)
524  zero=0
525  DO k=1, TOTAL_NUMBER_NODES_XI(3)
526    !Apex nodes
527    j=1
528    i=1
529    node_idx=node_idx+1
530    CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
531    IF(NodeDomain==ComputationalNodeNumber) THEN
532      FibreFieldAngle=(/zero,zero,zero/) 
533      DO component_idx=1,FieldFibreNumberOfComponents
534        CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
535          & DerivativeUserNumber,node_idx,component_idx,FibreFieldAngle(component_idx),Err)
536      ENDDO
537    ENDIF
538    theta=atan(FIBRE_SLOPE_CHANGE*XI3+FIBRE_SLOPE_INTERSECTION)
539    DO j=2, TOTAL_NUMBER_NODES_XI(2) 
540      nu=PI-XI2delta*(j-1)
541      omega=PI/2+cos(2*nu)*atan(SHEET_SLOPE_BASE_ENDO)*(-2*XI3+1)
542      DO i=1, TOTAL_NUMBER_NODES_XI(1)
543        node_idx=node_idx+1
544        CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
545        IF(NodeDomain==ComputationalNodeNumber) THEN
546          FibreFieldAngle=(/theta,zero,omega/)
547          DO component_idx=1,FieldFibreNumberOfComponents
548            CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
549              & DerivativeUserNumber, node_idx,component_idx,FibreFieldAngle(component_idx),Err)
550          ENDDO
551        ENDIF
552      ENDDO
553    ENDDO
554    XI3=XI3+XI3delta
555  ENDDO
556
557  !Create a material field and attach it to the geometric field
558  CALL CMISSField_Initialise(MaterialFieldSolid,Err)
559  CALL CMISSField_CreateStart(FieldMaterialUserNumber,Region,MaterialFieldSolid,Err)
560  CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err)
561  CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err)        
562  CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err)
563  CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialNumberOfVariables,Err)
564  CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialNumberOfComponents,Err)  
565  CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
566  CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
567  CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
568  CALL CMISSField_CreateFinish(MaterialFieldSolid,Err)
569
570  !Set Mooney-Rivlin constants 
571  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,C1,Err)
572  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,C2,Err)
573  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,C3,Err)
574
575  !
576  !================================================================================================================================
577  !
578
579  !EQUATIONS SETS
580  DO icompartment = 1,Ncompartments
581    EquationsSetFieldUserNumberDarcy = 100_CMISSIntg+icompartment
582    EquationsSetUserNumberDarcy = 200_CMISSIntg+icompartment
583  !Create the equations set for ALE Darcy
584    CALL CMISSField_Initialise(EquationsSetFieldDarcy(icompartment),Err)
585    CALL CMISSEquationsSet_Initialise(EquationsSetDarcy(icompartment),Err)
586    CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricFieldDarcy, &
587      & CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
588      & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
589      & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy(icompartment),EquationsSetDarcy(icompartment),Err)
590    !Finish creating the equations set
591    CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy(icompartment),Err)
592    !Set the values for the equations set field to be the current compartment number (1 - N), and the total number of compartments (N)
593    CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
594      & CMISS_FIELD_VALUES_SET_TYPE,1,icompartment,Err)
595    CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
596      & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
597  ENDDO
598
599  !--------------------------------------------------------------------------------------------------------------------------------
600  ! Solid
601
602  !Create the equations_set
603  CALL CMISSField_Initialise(EquationsSetFieldSolid,Err)
604  CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err)
605  CALL CMISSEquationsSet_CreateStart(EquationSetSolidUserNumber,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
606    & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
607    & EquationsSetFieldSolidUserNumber,EquationsSetFieldSolid,EquationsSetSolid,Err)
608  CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err)
609  !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)
610  !Need to store number of compartments, as finite elasticity uses this to calculate the total mass increase for the constiutive law
611  CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
612     & CMISS_FIELD_VALUES_SET_TYPE,1,0_CMISSIntg,Err)
613  CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
614     & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
615
616  !
617  !================================================================================================================================
618  !
619  ! Solid
620
621  !Create a dependent field with two variables and four components
622  CALL CMISSField_Initialise(DependentFieldSolid,Err)
623  !
624  CALL CMISSField_CreateStart(FieldDependentSolidUserNumber,Region,DependentFieldSolid,Err)
625  !
626  CALL CMISSField_TypeSet(DependentFieldSolid,CMISS_FIELD_GENERAL_TYPE,Err)
627  CALL CMISSField_MeshDecompositionSet(DependentFieldSolid,Decomposition,Err)
628  CALL CMISSField_GeometricFieldSet(DependentFieldSolid,GeometricFieldSolid,Err)
629  CALL CMISSField_DependentTypeSet(DependentFieldSolid,CMISS_FIELD_DEPENDENT_TYPE,Err)
630  !Create 2N+2 number of variables - 2 for solid, 2N for N Darcy compartments
631  FieldDependentSolidNumberOfVariables=2*Ncompartments+2
632  CALL CMISSField_NumberOfVariablesSet(DependentFieldSolid,FieldDependentSolidNumberOfVariables,Err)
633  !create two variables for each compartment
634  ALLOCATE(VariableTypes(2*Ncompartments+2))
635  DO num_var=1,Ncompartments+1
636     VariableTypes(2*num_var-1)=CMISS_FIELD_U_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
637     VariableTypes(2*num_var)=CMISS_FIELD_DELUDELN_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
638  ENDDO
639  CALL CMISSField_VariableTypesSet(DependentFieldSolid,VariableTypes,Err) 
640!   CALL CMISSField_VariableTypesSet(DependentFieldSolid,(/CMISS_FIELD_U_VARIABLE_TYPE, &
641!     & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err)
642    CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
643       & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
644    CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
645       & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
646  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
647  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
648    & FieldDependentSolidNumberOfComponents,Err)
649  DO icompartment=3,2*Ncompartments+2
650    CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
651  ENDDO
652!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
653!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
654!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
655
656  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidDisplMeshComponentNumber,Err)
657  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidDisplMeshComponentNumber,Err)
658  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidDisplMeshComponentNumber,Err)
659  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4, &
660    & CMISS_FIELD_NODE_BASED_INTERPOLATION, &
661    & Err)
662!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
663!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
664  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidLagrMultMeshComponentNumber,Err)
665!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
666!     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
667!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
668!     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
669!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
670!     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
671  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
672    & SolidDisplMeshComponentNumber, &
673    & Err)
674  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
675    & SolidDisplMeshComponentNumber, &
676    & Err)
677  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
678    & SolidDisplMeshComponentNumber, &
679    & Err)
680  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
681    & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
682!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
683!     & CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
684!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
685  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
686    & SolidLagrMultMeshComponentNumber, &
687    & Err)
688  !loop over the number of compartments
689  DO icompartment=3,2*Ncompartments+2
690!     CALL CMISSField_DimensionSet(DependentFieldSolid,VariableTypes(icompartment), &
691!        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
692    !CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
693    DO componentnum=1,FieldDependentFluidNumberOfComponents-1
694    !set dimension type
695!     CALL CMISSField_DimensionSet(DependentField,VariableTypes(icompartment), &
696!        & CMISS_FIELD_SCALAR_DIMENSION_TYPE,Err)
697      CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, &
698       & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
699      CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, & 
700         & DarcyVelMeshComponentNumber,Err)
701    ENDDO
702      CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment), &
703       & FieldDependentFluidNumberOfComponents, &
704       & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
705!     CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
706!       & FieldDependentFluidNumberOfComponents,MESH_COMPONENT_NUMBER_PRESSURE,Err)
707    CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
708      & FieldDependentFluidNumberOfComponents,DarcyMassIncreaseMeshComponentNumber,Err)
709    
710  ENDDO
711
712!   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
713!   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
714!   !For this equation type, MESH_COMPONENT_NUMBER_PRESSURE is actually the mass increase component as the pressure is taken from the solid equations
715!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
716!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
717!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
718!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
719!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
720!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
721!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
722!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
723
724  !
725  CALL CMISSField_ScalingTypeSet(DependentFieldSolid,CMISS_FIELD_UNIT_SCALING,Err)
726
727  CALL CMISSField_CreateFinish(DependentFieldSolid,Err)
728
729  !
730  !================================================================================================================================
731  !
732
733  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentSolidUserNumber,DependentFieldSolid,Err)
734  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err)
735
736  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialUserNumber,MaterialFieldSolid,Err)  
737  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err)
738
739  !
740  !================================================================================================================================
741  !
742  DO icompartment = 1,Ncompartments
743    CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy(icompartment),FieldDependentSolidUserNumber,&
744      & DependentFieldSolid,Err)
745    CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy(icompartment),Err)
746  ENDDO
747
748  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS+1
749    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
750      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
751    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
752      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
753    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U2_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
754      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
755    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U3_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
756      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
757  ENDDO
758  !
759  !================================================================================================================================
760  !
761  ALLOCATE(CouplingCoeffs(Ncompartments,Ncompartments))
762  IF(Ncompartments==2)THEN
763    CouplingCoeffs(1,1)=0.0E-01_CMISSDP
764!     CouplingCoeffs(1,2)=-1.0E-04_CMISSDP
765!     CouplingCoeffs(2,1)=-1.0E-04_CMISSDP
766    CouplingCoeffs(1,2)=0.0E-01_CMISSDP
767    CouplingCoeffs(2,1)=0.0E-01_CMISSDP
768    CouplingCoeffs(2,2)=0.0E-01_CMISSDP
769  ELSE IF(Ncompartments==3)THEN
770    CouplingCoeffs(1,1)=1.0E-02_CMISSDP
771    CouplingCoeffs(1,2)=1.0E-02_CMISSDP
772    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
773    CouplingCoeffs(2,1)=1.0E-02_CMISSDP
774    CouplingCoeffs(2,2)=2.0E-02_CMISSDP
775    CouplingCoeffs(2,3)=1.0E-02_CMISSDP
776    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
777    CouplingCoeffs(3,2)=1.0E-02_CMISSDP
778    CouplingCoeffs(3,3)=1.0E-02_CMISSDP
779  ELSE IF(Ncompartments==4)THEN
780    CouplingCoeffs(1,1)=0.0E-02_CMISSDP
781    CouplingCoeffs(1,2)=0.0E-02_CMISSDP
782    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
783    CouplingCoeffs(1,4)=0.0E-02_CMISSDP
784    CouplingCoeffs(2,1)=0.0E-02_CMISSDP
785    CouplingCoeffs(2,2)=0.0E-02_CMISSDP
786    CouplingCoeffs(2,3)=0.0E-02_CMISSDP
787    CouplingCoeffs(2,4)=0.0E-02_CMISSDP
788    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
789    CouplingCoeffs(3,2)=0.0E-02_CMISSDP
790    CouplingCoeffs(3,3)=0.0E-02_CMISSDP
791    CouplingCoeffs(3,4)=0.0E-02_CMISSDP
792    CouplingCoeffs(4,1)=0.0E-02_CMISSDP
793    CouplingCoeffs(4,2)=0.0E-02_CMISSDP
794    CouplingCoeffs(4,3)=0.0E-02_CMISSDP
795    CouplingCoeffs(4,4)=0.0E-02_CMISSDP
796  ELSE IF(Ncompartments==5)THEN
797    CouplingCoeffs(1,1)=0.0E-02_CMISSDP
798    CouplingCoeffs(1,2)=0.0E-02_CMISSDP
799    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
800    CouplingCoeffs(1,4)=0.0E-02_CMISSDP
801    CouplingCoeffs(1,5)=0.0E-02_CMISSDP
802    CouplingCoeffs(2,1)=0.0E-02_CMISSDP
803    CouplingCoeffs(2,2)=0.0E-02_CMISSDP
804    CouplingCoeffs(2,3)=0.0E-02_CMISSDP
805    CouplingCoeffs(2,4)=0.0E-02_CMISSDP
806    CouplingCoeffs(2,5)=0.0E-02_CMISSDP
807    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
808    CouplingCoeffs(3,2)=0.0E-02_CMISSDP
809    CouplingCoeffs(3,3)=0.0E-02_CMISSDP
810    CouplingCoeffs(3,4)=0.0E-02_CMISSDP
811    CouplingCoeffs(3,5)=0.0E-02_CMISSDP
812    CouplingCoeffs(4,1)=0.0E-02_CMISSDP
813    CouplingCoeffs(4,2)=0.0E-02_CMISSDP
814    CouplingCoeffs(4,3)=0.0E-02_CMISSDP
815    CouplingCoeffs(4,4)=0.0E-02_CMISSDP
816    CouplingCoeffs(4,5)=0.0E-02_CMISSDP
817    CouplingCoeffs(5,1)=0.0E-02_CMISSDP
818    CouplingCoeffs(5,2)=0.0E-02_CMISSDP
819    CouplingCoeffs(5,3)=0.0E-02_CMISSDP
820    CouplingCoeffs(5,4)=0.0E-02_CMISSDP
821    CouplingCoeffs(5,5)=0.0E-02_CMISSDP
822  ELSE
823    write(*,*) "Can't initialise coupling coefficients array."
824  ENDIF
825  !Define the material parameters for each compartments' constitutive law (for determining pressure)
826  Nparams=3
827  ALLOCATE(ConstitutiveParams(Ncompartments,Nparams))
828  IF(Ncompartments==2)THEN
829    ConstitutiveParams(1,1)=10.0E-01_CMISSDP
830    ConstitutiveParams(1,2)=10.0E-01_CMISSDP
831    ConstitutiveParams(1,3)=10.0E-01_CMISSDP
832    ConstitutiveParams(2,1)=10.0E-01_CMISSDP
833    ConstitutiveParams(2,2)=10.0E-01_CMISSDP
834    ConstitutiveParams(2,3)=10.0E-01_CMISSDP
835  ELSE IF(Ncompartments==3)THEN
836    ConstitutiveParams(1,1)=1.0E-02_CMISSDP
837    ConstitutiveParams(1,2)=1.0E-02_CMISSDP
838    ConstitutiveParams(1,3)=0.0E-02_CMISSDP
839    ConstitutiveParams(2,1)=1.0E-02_CMISSDP
840    ConstitutiveParams(2,2)=2.0E-02_CMISSDP
841    ConstitutiveParams(2,3)=1.0E-02_CMISSDP
842    ConstitutiveParams(3,1)=0.0E-02_CMISSDP
843    ConstitutiveParams(3,2)=1.0E-02_CMISSDP
844    ConstitutiveParams(3,3)=1.0E-02_CMISSDP
845  ELSE IF(Ncompartments==4)THEN
846    ConstitutiveParams(1,1)=0.0E-02_CMISSDP
847    ConstitutiveParams(1,2)=0.0E-02_CMISSDP
848    ConstitutiveParams(1,3)=0.0E-02_CMISSDP
849    ConstitutiveParams(2,1)=0.0E-02_CMISSDP
850    ConstitutiveParams(2,2)=0.0E-02_CMISSDP
851    ConstitutiveParams(2,3)=0.0E-02_CMISSDP
852    ConstitutiveParams(3,1)=0.0E-02_CMISSDP
853    ConstitutiveParams(3,2)=0.0E-02_CMISSDP
854    ConstitutiveParams(3,3)=0.0E-02_CMISSDP
855    ConstitutiveParams(4,1)=0.0E-02_CMISSDP
856    ConstitutiveParams(4,2)=0.0E-02_CMISSDP
857    ConstitutiveParams(4,3)=0.0E-02_CMISSDP
858  ELSE IF(Ncompartments==5)THEN
859    ConstitutiveParams(1,1)=0.0E-02_CMISSDP
860    ConstitutiveParams(1,2)=0.0E-02_CMISSDP
861    ConstitutiveParams(1,3)=0.0E-02_CMISSDP
862    ConstitutiveParams(2,1)=0.0E-02_CMISSDP
863    ConstitutiveParams(2,2)=0.0E-02_CMISSDP
864    ConstitutiveParams(2,3)=0.0E-02_CMISSDP
865    ConstitutiveParams(3,1)=0.0E-02_CMISSDP
866    ConstitutiveParams(3,2)=0.0E-02_CMISSDP
867    ConstitutiveParams(3,3)=0.0E-02_CMISSDP
868    ConstitutiveParams(4,1)=0.0E-02_CMISSDP
869    ConstitutiveParams(4,2)=0.0E-02_CMISSDP
870    ConstitutiveParams(4,3)=0.0E-02_CMISSDP
871    ConstitutiveParams(5,1)=0.0E-02_CMISSDP
872    ConstitutiveParams(5,2)=0.0E-02_CMISSDP
873    ConstitutiveParams(5,3)=0.0E-02_CMISSDP
874  ELSE
875    write(*,*) "Can't initialise constitutive parameters array."
876  ENDIF
877  !MATERIALS FIELDS
878  !Auto-created field contains a U variable type to store the diffusion coefficient(s)
879  !It also contains a V variable type to store the coupling coefficients 
880  DO icompartment = 1,Ncompartments
881    MaterialsFieldUserNumberDarcy = 400+icompartment
882    CALL CMISSField_Initialise(MaterialsFieldDarcy(icompartment),Err)
883    CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy(icompartment),MaterialsFieldUserNumberDarcy,&
884         & MaterialsFieldDarcy(icompartment),Err)
885    CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy(icompartment),Err)
886  END DO
887  DO icompartment = 1,Ncompartments
888    CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
889      & CMISS_FIELD_VALUES_SET_TYPE, &
890      & MaterialsFieldUserNumberDarcyPorosity,POROSITY_PARAM_DARCY,Err)
891    CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
892      & CMISS_FIELD_VALUES_SET_TYPE, &
893      & MaterialsFieldUserNumberDarcyPermOverVis,PERM_OVER_VIS_PARAM_DARCY,Err)
894  END DO
895  DO icompartment = 1, Ncompartments
896    DO COMPONENT_NUMBER=1, Ncompartments
897      CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
898         & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
899!         CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
900!           & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
901    END DO
902  END DO
903  DO icompartment = 1, Ncompartments
904    DO COMPONENT_NUMBER=1,Nparams
905      CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U1_VARIABLE_TYPE, &
906         & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,ConstitutiveParams(icompartment,COMPONENT_NUMBER),Err)
907!         CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
908!           & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
909    END DO
910  END DO
911
912  !
913  !================================================================================================================================
914  !
915
916  !EQUATIONS SET EQUATIONS
917
918  !Darcy
919  DO icompartment=1,Ncompartments
920    !Create the equations set equations
921    CALL CMISSEquations_Initialise(EquationsDarcy(icompartment),Err)
922    CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy(icompartment),EquationsDarcy(icompartment),Err)
923    !Set the equations matrices sparsity type
924    CALL CMISSEquations_SparsityTypeSet(EquationsDarcy(icompartment

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