/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 are truncated click here to view the full file

    

  1. !> \file
    2. 

  2. !> \author Christian Michler
    3. 

  3. !> \brief This is an example program to solve a coupled Finite Elastiticity Darcy equation on a cylindrical geometry using openCMISS calls.
    4. 

  4. !>
    5. 

  5. !> \section LICENSE
    6. 

  6. !>
    7. 

  7. !> Version: MPL 1.1/GPL 2.0/LGPL 2.1
    8. 

  8. !>
    9. 

  9. !> The contents of this file are subject to the Mozilla Public License
    10. 

  10. !> Version 1.1 (the "License"); you may not use this file except in
    11. 

  11. !> compliance with the License. You may obtain a copy of the License at
    12. 

  12. !> http://www.mozilla.org/MPL/
    13. 

  13. !>
    14. 

  14. !> Software distributed under the License is distributed on an "AS IS"
    15. 

  15. !> basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
    16. 

  16. !> License for the specific language governing rights and limitations
    17. 

  17. !> under the License.
    18. 

  18. !>
    19. 

  19. !> The Original Code is OpenCMISS
    20. 

  20. !>
    21. 

  21. !> The Initial Developer of the Original Code is University of Auckland,
    22. 

  22. !> Auckland, New Zealand and University of Oxford, Oxford, United
    23. 

  23. !> Kingdom. Portions created by the University of Auckland and University
    24. 

  24. !> of Oxford are Copyright (C) 2007 by the University of Auckland and
    25. 

  25. !> the University of Oxford. All Rights Reserved.
    26. 

  26. !>
    27. 

  27. !> Contributor(s): Christian Michler, Jack Lee
    28. 

  28. !>
    29. 

  29. !> Alternatively, the contents of this file may be used under the terms of
    30. 

  30. !> either the GNU General Public License Version 2 or later (the "GPL"), or
    31. 

  31. !> the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
    32. 

  32. !> in which case the provisions of the GPL or the LGPL are applicable instead
    33. 

  33. !> of those above. If you wish to allow use of your version of this file only
    34. 

  34. !> under the terms of either the GPL or the LGPL, and not to allow others to
    35. 

  35. !> use your version of this file under the terms of the MPL, indicate your
    36. 

  36. !> decision by deleting the provisions above and replace them with the notice
    37. 

  37. !> and other provisions required by the GPL or the LGPL. If you do not delete
    38. 

  38. !> the provisions above, a recipient may use your version of this file under
    39. 

  39. !> the terms of any one of the MPL, the GPL or the LGPL.
    40. 

  40. !>
    41. 

  41. 

  42. !> \example MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenDarcy/src/QuadraticEllipsoidDrivenDarcyExample.f90
    43. 

  43. !! Example program to solve coupled FiniteElasticityDarcy equations using OpenCMISS calls.
    44. 

  44. !! \par Latest Builds:
    45. 

  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. 

  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. 

  47. !!
    48. 

  48. !<
    49. 

  49. 

  50. ! !
    51. 

  51. ! !  This example considers a coupled Finite Elasticity Darcy problem on a cylindrical geometry
    52. 

  52. ! !
    53. 

  53. 

  54. !> Main program
    55. 

  55. 

  56. PROGRAM QUADRATICELLIPSOIDDRIVENMULTICOMPDARCYEXAMPLE
    57. 

  57. 

  58.   !
    59. 

  59.   !================================================================================================================================
    60. 

  60.   !
    61. 

  61. 

  62.   !PROGRAM LIBRARIES
    63. 

  63. 

  64.   USE OPENCMISS
    65. 

  65.   USE MPI
    66. 

  66. 

  67. #ifdef WIN32
    68. 

  68.   USE IFQWINCMISS
    69. 

  69. #endif
    70. 

  70. 

  71.   !
    72. 

  72.   !================================================================================================================================
    73. 

  73.   !
    74. 

  74. 

  75.   !PROGRAM VARIABLES AND TYPES
    76. 

  76. 

  77.   IMPLICIT NONE
    78. 

  78. 

  79.   !Test program parameters
    80. 

  80. 

  81.   REAL(CMISSDP), PARAMETER :: PI=3.14159_CMISSDP
    82. 

  82.   REAL(CMISSDP), PARAMETER :: LONG_AXIS=2.0_CMISSDP
    83. 

  83.   REAL(CMISSDP), PARAMETER :: SHORT_AXIS=1.0_CMISSDP
    84. 

  84.   REAL(CMISSDP), PARAMETER :: WALL_THICKNESS=0.5_CMISSDP
    85. 

  85.   REAL(CMISSDP), PARAMETER :: CUTOFF_ANGLE=1.5708_CMISSDP
    86. 

  86.   REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_INTERSECTION=1.73205_CMISSDP !Slope of fibres in base endocardium = 60 degrees
    87. 

  87.   REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_CHANGE=-3.4641_CMISSDP !Slope change of fibres from 60 to -60 degrees in transmural direction 
    88. 

  88.   REAL(CMISSDP), PARAMETER :: SHEET_SLOPE_BASE_ENDO=1.0_CMISSDP !Slope of sheet at base endocardium 
    89. 

  89. 

  90.   REAL(CMISSDP), PARAMETER :: INNER_PRESSURE=2.0_CMISSDP  !Positive is compressive
    91. 

  91.   REAL(CMISSDP), PARAMETER :: OUTER_PRESSURE=0.0_CMISSDP !Positive is compressive
    92. 

  92.   REAL(CMISSDP), PARAMETER :: C1= 2.0_CMISSDP
    93. 

  93.   REAL(CMISSDP), PARAMETER :: C2= 6.0_CMISSDP
    94. 

  94.   REAL(CMISSDP), PARAMETER :: C3=10.0_CMISSDP
    95. 

  95. 

  96.   INTEGER(CMISSIntg), PARAMETER :: NumberGlobalXElements=4  ! X ==NUMBER_GLOBAL_CIRCUMFERENTIAL_ELEMENTS
    97. 

  97.   INTEGER(CMISSIntg), PARAMETER :: NumberGlobalYElements=4  ! Y ==NUMBER_GLOBAL_LONGITUDINAL_ELEMENTS
    98. 

  98.   INTEGER(CMISSIntg), PARAMETER :: NumberGlobalZElements=1  ! Z ==NUMBER_GLOBAL_TRANSMURAL_ELEMENTS
    99. 

  99.   INTEGER(CMISSIntg) :: NumberOfDomains
    100. 

  100. 

  101.   INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
    102. 

  102.   INTEGER(CMISSIntg), PARAMETER :: NumberOfSpatialCoordinates=3
    103. 

  103.   INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=1
    104. 

  104.   INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=1
    105. 

  105.   INTEGER(CMISSIntg), PARAMETER :: QuadraticCollapsedBasisUserNumber=2
    106. 

  106.   INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=3
    107. 

  107.   INTEGER(CMISSIntg), PARAMETER :: LinearCollapsedBasisUserNumber=4
    108. 

  108.   INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=1
    109. 

  109.   INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=2
    110. 

  110.   INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=1
    111. 

  111.   INTEGER(CMISSIntg), PARAMETER :: DerivativeUserNumber=1
    112. 

  112. 

  113.   INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshDimensions=3
    114. 

  114.   INTEGER(CMISSIntg), PARAMETER :: NumberOfXiCoordinates=3
    115. 

  115.   INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshComponents=2
    116. 

  116.   INTEGER(CMISSIntg), PARAMETER :: QuadraticMeshComponentNumber=1
    117. 

  117.   INTEGER(CMISSIntg), PARAMETER :: LinearMeshComponentNumber=2
    118. 

  118.   INTEGER(CMISSIntg), PARAMETER :: TotalNumberOfElements=1
    119. 

  119. 

  120.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumberSolid=1
    121. 

  121.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumberDarcy=11
    122. 

  122.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1
    123. 

  123.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
    124. 

  124. 

  125.   INTEGER(CMISSIntg), PARAMETER :: FieldFibreUserNumber=2
    126. 

  126.   INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfVariables=1
    127. 

  127.   INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfComponents=3
    128. 

  128. 

  129.   INTEGER(CMISSIntg), PARAMETER :: FieldMaterialUserNumber=3
    130. 

  130.   INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfVariables=1
    131. 

  131.   INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfComponents=3 !2
    132. 

  132. 

  133.   INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidUserNumber=4
    134. 

  134.   INTEGER(CMISSIntg) :: FieldDependentSolidNumberOfVariables
    135. 

  135.   INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4
    136. 

  136.   INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4
    137. 

  137. 

  138.   INTEGER(CMISSIntg), PARAMETER :: EquationSetSolidUserNumber=55
    139. 

  139.   INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldSolidUserNumber=25
    140. 

  140.   INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=1
    141. 

  141.   INTEGER(CMISSIntg) :: EquationsSetFieldUserNumberDarcy
    142. 

  142.   INTEGER(CMISSIntg) :: icompartment,Ncompartments,num_var,componentnum,Nparams
    143. 

  143. 

  144.   INTEGER(CMISSIntg) :: MaterialsFieldUserNumberDarcy
    145. 

  145.   INTEGER(CMISSIntg) :: EquationsSetUserNumberDarcy
    146. 

  146. 

  147. 

  148.   INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1
    149. 

  149.   INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2
    150. 

  150.   INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1
    151. 

  151.   INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1
    152. 

  152.   INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1
    153. 

  153.   INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPorosity=1
    154. 

  154.   INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPermOverVis=2
    155. 

  155.   INTEGER(CMISSIntg), PARAMETER :: SolidDisplMeshComponentNumber=1
    156. 

  156.   INTEGER(CMISSIntg), PARAMETER :: SolidLagrMultMeshComponentNumber=2
    157. 

  157.   INTEGER(CMISSIntg), PARAMETER :: SolidGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
    158. 

  158.   !Program types
    159. 

  159. 

  160. 

  161.   !Program variables
    162. 

  162. 

  163.   INTEGER(CMISSIntg) :: MPI_IERROR
    164. 

  164.   INTEGER(CMISSIntg) :: EquationsSetIndex  
    165. 

  165.   INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber
    166. 

  166.   REAL(CMISSDP) :: FibreFieldAngle(3) 
    167. 

  167.   REAL(CMISSDP) :: nu,theta,omega,XI3,XI3delta,XI2delta, zero
    168. 

  168.   INTEGER(CMISSIntg) ::i,j,k,component_idx,node_idx,TOTAL_NUMBER_NODES_XI(3)
    169. 

  169.   !For grabbing surfaces
    170. 

  170.   INTEGER(CMISSIntg) :: InnerNormalXi,OuterNormalXi,TopNormalXi
    171. 

  171.   INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodes(:)
    172. 

  172.   INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodes(:)
    173. 

  173.   INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodes(:)
    174. 

  174.   INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodesDarcyVel(:)
    175. 

  175.   INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodesDarcyVel(:)
    176. 

  176.   INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodesDarcyVel(:)
    177. 

  177.   INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
    178. 

  178.   REAL(CMISSDP) :: XCoord,YCoord,ZCoord
    179. 

  179.   LOGICAL :: X_FIXED,Y_FIXED,X_OKAY,Y_OKAY
    180. 

  180. 

  181.   INTEGER(CMISSIntg) :: GeometricFieldDarcyMeshComponentNumber, DarcyVelMeshComponentNumber, DarcyMassIncreaseMeshComponentNumber
    182. 

  182.   INTEGER(CMISSIntg) :: MeshComponentNumber_dummy
    183. 

  183. 

  184.   INTEGER(CMISSIntg) :: NUMBER_OF_DOMAINS
    185. 

  185. 

  186.   INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
    187. 

  187. 

  188.   INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
    189. 

  189.   INTEGER(CMISSIntg) :: RESTART_VALUE
    190. 

  190. 

  191.   INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT
    192. 

  192.   INTEGER(CMISSIntg) :: COMPONENT_NUMBER
    193. 

  193.   INTEGER(CMISSIntg) :: CONDITION
    194. 

  194. 

  195.   INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY
    196. 

  196.   INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE
    197. 

  197.   INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE
    198. 

  198. 

  199.   REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z
    200. 

  200.   REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2
    201. 

  201.   REAL(CMISSDP) :: GEOMETRY_TOLERANCE
    202. 

  202.   INTEGER(CMISSIntg) :: EDGE_COUNT
    203. 

  203.   INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID
    204. 

  204.   REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4)
    205. 

  205.   REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
    206. 

  206.   REAL(CMISSDP) :: RELATIVE_TOLERANCE
    207. 

  207.   REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
    208. 

  208.   REAL(CMISSDP) :: LINESEARCH_ALPHA
    209. 

  209.   REAL(CMISSDP) :: VALUE
    210. 

  210.   REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY
    211. 

  211. 

  212.   LOGICAL :: EXPORT_FIELD_IO
    213. 

  213.   LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG
    214. 

  214. 

  215.   !CMISS variables
    216. 

  216. 

  217.   TYPE(CMISSBasisType) :: QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis
    218. 

  218.   TYPE(CMISSBoundaryConditionsType) :: BoundaryConditions
    219. 

  219.   TYPE(CMISSCoordinateSystemType) :: CoordinateSystem, WorldCoordinateSystem
    220. 

  220.   TYPE(CMISSMeshType) :: Mesh
    221. 

  221.   TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
    222. 

  222.   TYPE(CMISSDecompositionType) :: Decomposition
    223. 

  223.   TYPE(CMISSEquationsType) :: Equations
    224. 

  224.   TYPE(CMISSEquationsSetType) :: EquationsSetSolid
    225. 

  225.   TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid
    226. 

  226.   TYPE(CMISSFieldType) :: DependentFieldSolid,EquationsSetFieldSolid
    227. 

  227. 

  228.   TYPE(CMISSFieldType) :: GeometricFieldDarcy
    229. 

  229. 

  230.   TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: MaterialsFieldDarcy
    231. 

  231.   TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: EquationsSetFieldDarcy
    232. 

  232.   !Boundary conditions
    233. 

  233.   TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy
    234. 

  234.   !Equations sets
    235. 

  235.   TYPE(CMISSEquationsSetType), ALLOCATABLE, DIMENSION(:) :: EquationsSetDarcy
    236. 

  236.   !Equations
    237. 

  237.   TYPE(CMISSEquationsType), ALLOCATABLE, DIMENSION(:) :: EquationsDarcy
    238. 

  238. 

  239. 

  240. 

  241.   TYPE(CMISSFieldsType) :: Fields
    242. 

  242.   TYPE(CMISSProblemType) :: Problem
    243. 

  243.   TYPE(CMISSRegionType) :: Region,WorldRegion
    244. 

  244.   TYPE(CMISSSolverType) :: SolverSolid,LinearSolverSolid
    245. 

  245.   TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid
    246. 

  246.   TYPE(CMISSControlLoopType) :: ControlLoop
    247. 

  247. 

  248.   !Solvers
    249. 

  249.   TYPE(CMISSSolverType) :: DynamicSolverDarcy
    250. 

  250.   TYPE(CMISSSolverType) :: LinearSolverDarcy
    251. 

  251.   !Solver equations
    252. 

  252.   TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy
    253. 

  253. 

  254.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME
    255. 

  255.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME
    256. 

  256.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA
    257. 

  257.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
    258. 

  258. 

  259. #ifdef WIN32
    260. 

  260.   !Quickwin type
    261. 

  261.   LOGICAL :: QUICKWIN_STATUS=.FALSE.
    262. 

  262.   TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
    263. 

  263. #endif
    264. 

  264. 

  265.   !Generic CMISS variables
    266. 

  266.   INTEGER(CMISSIntg) :: Err
    267. 

  267.   INTEGER(CMISSIntg), ALLOCATABLE, DIMENSION(:) :: VariableTypes
    268. 

  268.   REAL(CMISSDP), ALLOCATABLE, DIMENSION(:,:) :: CouplingCoeffs,ConstitutiveParams
    269. 

  269. #ifdef WIN32
    270. 

  270.   !Initialise QuickWin
    271. 

  271.   QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
    272. 

  272.   QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
    273. 

  273.   QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
    274. 

  274.   !Set the window parameters
    275. 

  275.   QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
    276. 

  276.   !If attempt fails set with system estimated values
    277. 

  277.   IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
    278. 

  278. #endif
    279. 

  279. 

  280.   !Set initial values
    281. 

  281.   NUMBER_OF_DIMENSIONS=3_CMISSIntg
    282. 

  282.   INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
    283. 

  283.   INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
    284. 

  284.   INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
    285. 

  285.   INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
    286. 

  286.   !Set material parameters
    287. 

  287.   POROSITY_PARAM_DARCY=0.1_CMISSDP
    288. 

  288.   PERM_OVER_VIS_PARAM_DARCY=1.0_CMISSDP
    289. 

  289.   !Set output parameter
    290. 

  290.   !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
    291. 

  291.   DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
    292. 

  292.   LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_SOLVER_OUTPUT
    293. 

  293.   !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
    294. 

  294.   EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
    295. 

  295. 

  296.   !Set time parameter
    297. 

  297.   DYNAMIC_SOLVER_DARCY_START_TIME=0.0_CMISSDP
    298. 

  298.   DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP
    299. 

  299.   DYNAMIC_SOLVER_DARCY_STOP_TIME=5_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
    300. 

  300.   DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
    301. 

  301.   !Set result output parameter
    302. 

  302.   DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
    303. 

  303.   !Set solver parameters
    304. 

  304.   LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
    305. 

  305.   RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
    306. 

  306.   ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
    307. 

  307.   DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
    308. 

  308.   MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
    309. 

  309.   RESTART_VALUE=30_CMISSIntg !default: 30
    310. 

  310.   LINESEARCH_ALPHA=1.0_CMISSDP
    311. 

  311.   icompartment =1_CMISSIntg
    312. 

  312.   Ncompartments=4_CMISSIntg
    313. 

  313. 

  314.   !LinearMeshComponentNumber/QuadraticMeshComponentNumber
    315. 

  315.   DarcyVelMeshComponentNumber = LinearMeshComponentNumber
    316. 

  316.   DarcyMassIncreaseMeshComponentNumber = LinearMeshComponentNumber
    317. 

  317. !   GeometricFieldDarcyMeshComponentNumber = DarcyVelMeshComponentNumber
    318. 

  318.   GeometricFieldDarcyMeshComponentNumber = QuadraticMeshComponentNumber
    319. 

  319. 

  320.   ALLOCATE (EquationsSetDarcy(Ncompartments))
    321. 

  321.   ALLOCATE (EquationsSetFieldDarcy(Ncompartments))
    322. 

  322.   ALLOCATE (MaterialsFieldDarcy(Ncompartments))
    323. 

  323.   ALLOCATE (EquationsDarcy(Ncompartments))
    324. 

  324.   !
    325. 

  325.   !================================================================================================================================
    326. 

  326.   !
    327. 

  327. 

  328.   !Intialise cmiss
    329. 

  329.   CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
    330. 

  330. 

  331.   CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
    332. 

  332. 

  333.   WRITE(*,'(A)') "Program starting."
    334. 

  334. 

  335.   !Set all diganostic levels on for testing
    336. 

  336.   CALL CMISSDiagnosticsSetOn(CMISS_FROM_DIAG_TYPE,(/1,2,3,4,5/),"Diagnostics",(/"PROBLEM_FINITE_ELEMENT_CALCULATE"/),Err)
    337. 

  337. 

  338.   !Get the number of computational nodes and this computational node number
    339. 

  339.   CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
    340. 

  340.   CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
    341. 

  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. 

  343. !   CALL MPI_BCAST(NumberGlobalXElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
    344. 

  344. !   CALL MPI_BCAST(NumberGlobalYElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
    345. 

  345. !   CALL MPI_BCAST(NumberGlobalZElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
    346. 

  346. !   CALL MPI_BCAST(NumberOfDomains,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
    347. 

  347.   NumberOfDomains=NumberOfComputationalNodes
    348. 

  348. 

  349.   !Create a CS - default is 3D rectangular cartesian CS with 0,0,0 as origin
    350. 

  350.   CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
    351. 

  351.   CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
    352. 

  352.   CALL CMISSCoordinateSystem_TypeSet(CoordinateSystem,CMISS_COORDINATE_RECTANGULAR_CARTESIAN_TYPE,Err)
    353. 

  353.   CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NumberOfSpatialCoordinates,Err)
    354. 

  354.   CALL CMISSCoordinateSystem_OriginSet(CoordinateSystem,(/0.0_CMISSDP,0.0_CMISSDP,0.0_CMISSDP/),Err)
    355. 

  355.   CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
    356. 

  356. 

  357.   !
    358. 

  358.   !================================================================================================================================
    359. 

  359.   !
    360. 

  360. 

  361.   !Create a region and assign the CS to the region
    362. 

  362.   CALL CMISSRegion_Initialise(Region,Err)
    363. 

  363.   CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
    364. 

  364.   CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
    365. 

  365.   CALL CMISSRegion_CreateFinish(Region,Err)
    366. 

  366. 

  367.   !
    368. 

  368.   !================================================================================================================================
    369. 

  369.   !
    370. 

  370. 

  371.   !Define basis functions - tri-linear Lagrange and tri-Quadratic Lagrange, each with collapsed variant
    372. 

  372.     !Quadratic Basis
    373. 

  373.   CALL CMISSBasis_Initialise(QuadraticBasis,Err)
    374. 

  374.   CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
    375. 

  375.   CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
    376. 

  376.     & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
    377. 

  377.   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
    378. 

  378.     & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
    379. 

  379. !     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
    380. 

  380.   CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err) !Have to do this
    381. 

  381.   CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
    382. 

  382. 

  383.     !Collapsed Quadratic Basis
    384. 

  384.   CALL CMISSBasis_Initialise(QuadraticCollapsedBasis,Err)
    385. 

  385.   CALL CMISSBasis_CreateStart(QuadraticCollapsedBasisUserNumber,QuadraticCollapsedBasis,Err)
    386. 

  386.   CALL CMISSBasis_TypeSet(QuadraticCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
    387. 

  387.   CALL CMISSBasis_NumberOfXiSet(QuadraticCollapsedBasis,NumberOfXiCoordinates,Err)
    388. 

  388.   CALL CMISSBasis_InterpolationXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
    389. 

  389.        & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
    390. 

  390.   CALL CMISSBasis_CollapsedXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED, &
    391. 

  391.        & CMISS_BASIS_COLLAPSED_AT_XI0,CMISS_BASIS_NOT_COLLAPSED/),Err)
    392. 

  392.   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticCollapsedBasis, &
    393. 

  393.        & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)  
    394. 

  394. !     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
    395. 

  395.   CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticCollapsedBasis,.true.,Err) !Have to do this
    396. 

  396.   CALL CMISSBasis_CreateFinish(QuadraticCollapsedBasis,Err)
    397. 

  397. 

  398.     !Linear Basis
    399. 

  399.   CALL CMISSBasis_Initialise(LinearBasis,Err)
    400. 

  400.   CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
    401. 

  401.   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
    402. 

  402.     & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
    403. 

  403. !     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
    404. 

  404.   CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
    405. 

  405.   CALL CMISSBasis_CreateFinish(LinearBasis,Err)
    406. 

  406. 

  407.     !Collapsed Linear Basis
    408. 

  408.   CALL CMISSBasis_Initialise(LinearCollapsedBasis,Err)
    409. 

  409.   CALL CMISSBasis_CreateStart(LinearCollapsedBasisUserNumber,LinearCollapsedBasis,Err)
    410. 

  410.   CALL CMISSBasis_TypeSet(LinearCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
    411. 

  411.   CALL CMISSBasis_NumberOfXiSet(LinearCollapsedBasis,NumberOfXiCoordinates,Err)
    412. 

  412.   CALL CMISSBasis_InterpolationXiSet(LinearCollapsedBasis,(/CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION, &
    413. 

  413.        & CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION,CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION/),Err)
    414. 

  414.   CALL CMISSBasis_CollapsedXiSet(LinearCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED,CMISS_BASIS_COLLAPSED_AT_XI0, &
    415. 

  415.     & CMISS_BASIS_NOT_COLLAPSED/),Err)
    416. 

  416.   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearCollapsedBasis, &
    417. 

  417.        & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
    418. 

  418. !     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
    419. 

  419.   CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearCollapsedBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
    420. 

  420.   CALL CMISSBasis_CreateFinish(LinearCollapsedBasis,Err)
    421. 

  421. 

  422.   !
    423. 

  423.   !================================================================================================================================
    424. 

  424.   !
    425. 

  425. 

  426.   !Start the creation of a generated ellipsoid mesh
    427. 

  427.   CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
    428. 

  428.   CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
    429. 

  429.   !Set up an ellipsoid mesh
    430. 

  430.   CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_MESH_TYPE,Err)
    431. 

  431.   !Set the quadratic and linear bases
    432. 

  432.   CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,[QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis],Err)
    433. 

  433.   !Define the mesh on the region
    434. 

  434.   CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/LONG_AXIS,SHORT_AXIS,WALL_THICKNESS,CUTOFF_ANGLE/),Err)
    435. 

  435.   CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NumberGlobalXElements,NumberGlobalYElements, &
    436. 

  436.     & NumberGlobalZElements/),Err)
    437. 

  437.   
    438. 

  438.   !Finish the creation of a generated mesh in the region
    439. 

  439.   CALL CMISSMesh_Initialise(Mesh,Err)
    440. 

  440.   CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
    441. 

  441. 

  442.   !
    443. 

  443.   !================================================================================================================================
    444. 

  444.   !
    445. 

  445. 

  446.   !Create a decomposition
    447. 

  447.   CALL CMISSDecomposition_Initialise(Decomposition,Err)
    448. 

  448.   CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
    449. 

  449.   CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
    450. 

  450.   CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
    451. 

  451.   CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.TRUE.,Err)
    452. 

  452.   CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
    453. 

  453. 

  454.   !
    455. 

  455.   !================================================================================================================================
    456. 

  456.   !
    457. 

  457. 

  458.   ! --- GeometricFieldSolid ---
    459. 

  459.   !Create a field to put the geometry (default is geometry)
    460. 

  460.   CALL CMISSField_Initialise(GeometricFieldSolid,Err)
    461. 

  461.   CALL CMISSField_CreateStart(FieldGeometryUserNumberSolid,Region,GeometricFieldSolid,Err)
    462. 

  462.   CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err)
    463. 

  463.   CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err)
    464. 

  464.   CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometryNumberOfVariables,Err)
    465. 

  465.   CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)  
    466. 

  466.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
    467. 

  467.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
    468. 

  468.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
    469. 

  469.   CALL CMISSField_CreateFinish(GeometricFieldSolid,Err)
    470. 

  470. 

  471.   !Update the geometric field parameters
    472. 

  472.   CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err)
    473. 

  473. 

  474.   !
    475. 

  475.   !================================================================================================================================
    476. 

  476.   !
    477. 

  477. 

  478.   ! --- GeometricFieldDarcy ---
    479. 

  479.   !Create a field to put the geometry (default is geometry)
    480. 

  480.   CALL CMISSField_Initialise(GeometricFieldDarcy,Err)
    481. 

  481.   CALL CMISSField_CreateStart(FieldGeometryUserNumberDarcy,Region,GeometricFieldDarcy,Err)
    482. 

  482.   CALL CMISSField_MeshDecompositionSet(GeometricFieldDarcy,Decomposition,Err)
    483. 

  483.   CALL CMISSField_TypeSet(GeometricFieldDarcy,CMISS_FIELD_GEOMETRIC_TYPE,Err)
    484. 

  484.   CALL CMISSField_NumberOfVariablesSet(GeometricFieldDarcy,FieldGeometryNumberOfVariables,Err)
    485. 

  485.   CALL CMISSField_NumberOfComponentsSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)  
    486. 

  486.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,1, &
    487. 

  487.     & GeometricFieldDarcyMeshComponentNumber,Err)
    488. 

  488.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,2, &
    489. 

  489.     & GeometricFieldDarcyMeshComponentNumber,Err)
    490. 

  490.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,3, &
    491. 

  491.     & GeometricFieldDarcyMeshComponentNumber,Err)
    492. 

  492.   CALL CMISSField_CreateFinish(GeometricFieldDarcy,Err)
    493. 

  493. 

  494.   !Update the geometric field parameters
    495. 

  495.   CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldDarcy,Err)
    496. 

  496. 

  497.   !
    498. 

  498.   !================================================================================================================================
    499. 

  499.   !
    500. 

  500. 

  501.   !Create a fibre field and attach it to the geometric field
    502. 

  502.   CALL CMISSField_Initialise(FibreFieldSolid,Err)
    503. 

  503.   CALL CMISSField_CreateStart(FieldFibreUserNumber,Region,FibreFieldSolid,Err)
    504. 

  504.   CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err)
    505. 

  505.   CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err)        
    506. 

  506.   CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err)
    507. 

  507.   CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreNumberOfVariables,Err)
    508. 

  508.   CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreNumberOfComponents,Err)  
    509. 

  509.   CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
    510. 

  510.   CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
    511. 

  511.   CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
    512. 

  512.   CALL CMISSField_CreateFinish(FibreFieldSolid,Err)
    513. 

  513. 

  514.   !Set Fibre directions (this block is parallel-untested)
    515. 

  515.   node_idx=0  
    516. 

  516.   !This is valid only for quadratic basis functions
    517. 

  517.   TOTAL_NUMBER_NODES_XI(1)=NumberGlobalXElements*2
    518. 

  518.   TOTAL_NUMBER_NODES_XI(2)=NumberGlobalYElements*2+1
    519. 

  519.   TOTAL_NUMBER_NODES_XI(3)=NumberGlobalZElements*2+1
    520. 

  520. 

  521.   XI2delta=(PI-CUTOFF_ANGLE)/(TOTAL_NUMBER_NODES_XI(2)-1)
    522. 

  522.   XI3=0
    523. 

  523.   XI3delta=(1.0)/(TOTAL_NUMBER_NODES_XI(3)-1)
    524. 

  524.   zero=0
    525. 

  525.   DO k=1, TOTAL_NUMBER_NODES_XI(3)
    526. 

  526.     !Apex nodes
    527. 

  527.     j=1
    528. 

  528.     i=1
    529. 

  529.     node_idx=node_idx+1
    530. 

  530.     CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
    531. 

  531.     IF(NodeDomain==ComputationalNodeNumber) THEN
    532. 

  532.       FibreFieldAngle=(/zero,zero,zero/) 
    533. 

  533.       DO component_idx=1,FieldFibreNumberOfComponents
    534. 

  534.         CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
    535. 

  535.           & DerivativeUserNumber,node_idx,component_idx,FibreFieldAngle(component_idx),Err)
    536. 

  536.       ENDDO
    537. 

  537.     ENDIF
    538. 

  538.     theta=atan(FIBRE_SLOPE_CHANGE*XI3+FIBRE_SLOPE_INTERSECTION)
    539. 

  539.     DO j=2, TOTAL_NUMBER_NODES_XI(2) 
    540. 

  540.       nu=PI-XI2delta*(j-1)
    541. 

  541.       omega=PI/2+cos(2*nu)*atan(SHEET_SLOPE_BASE_ENDO)*(-2*XI3+1)
    542. 

  542.       DO i=1, TOTAL_NUMBER_NODES_XI(1)
    543. 

  543.         node_idx=node_idx+1
    544. 

  544.         CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err)
    545. 

  545.         IF(NodeDomain==ComputationalNodeNumber) THEN
    546. 

  546.           FibreFieldAngle=(/theta,zero,omega/)
    547. 

  547.           DO component_idx=1,FieldFibreNumberOfComponents
    548. 

  548.             CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
    549. 

  549.               & DerivativeUserNumber, node_idx,component_idx,FibreFieldAngle(component_idx),Err)
    550. 

  550.           ENDDO
    551. 

  551.         ENDIF
    552. 

  552.       ENDDO
    553. 

  553.     ENDDO
    554. 

  554.     XI3=XI3+XI3delta
    555. 

  555.   ENDDO
    556. 

  556. 

  557.   !Create a material field and attach it to the geometric field
    558. 

  558.   CALL CMISSField_Initialise(MaterialFieldSolid,Err)
    559. 

  559.   CALL CMISSField_CreateStart(FieldMaterialUserNumber,Region,MaterialFieldSolid,Err)
    560. 

  560.   CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err)
    561. 

  561.   CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err)        
    562. 

  562.   CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err)
    563. 

  563.   CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialNumberOfVariables,Err)
    564. 

  564.   CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialNumberOfComponents,Err)  
    565. 

  565.   CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
    566. 

  566.   CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
    567. 

  567.   CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
    568. 

  568.   CALL CMISSField_CreateFinish(MaterialFieldSolid,Err)
    569. 

  569. 

  570.   !Set Mooney-Rivlin constants 
    571. 

  571.   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,C1,Err)
    572. 

  572.   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,C2,Err)
    573. 

  573.   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,C3,Err)
    574. 

  574. 

  575.   !
    576. 

  576.   !================================================================================================================================
    577. 

  577.   !
    578. 

  578. 

  579.   !EQUATIONS SETS
    580. 

  580.   DO icompartment = 1,Ncompartments
    581. 

  581.     EquationsSetFieldUserNumberDarcy = 100_CMISSIntg+icompartment
    582. 

  582.     EquationsSetUserNumberDarcy = 200_CMISSIntg+icompartment
    583. 

  583.   !Create the equations set for ALE Darcy
    584. 

  584.     CALL CMISSField_Initialise(EquationsSetFieldDarcy(icompartment),Err)
    585. 

  585.     CALL CMISSEquationsSet_Initialise(EquationsSetDarcy(icompartment),Err)
    586. 

  586.     CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricFieldDarcy, &
    587. 

  587.       & CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
    588. 

  588.       & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
    589. 

  589.       & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy(icompartment),EquationsSetDarcy(icompartment),Err)
    590. 

  590.     !Finish creating the equations set
    591. 

  591.     CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy(icompartment),Err)
    592. 

  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. 

  593.     CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    594. 

  594.       & CMISS_FIELD_VALUES_SET_TYPE,1,icompartment,Err)
    595. 

  595.     CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    596. 

  596.       & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
    597. 

  597.   ENDDO
    598. 

  598. 

  599.   !--------------------------------------------------------------------------------------------------------------------------------
    600. 

  600.   ! Solid
    601. 

  601. 

  602.   !Create the equations_set
    603. 

  603.   CALL CMISSField_Initialise(EquationsSetFieldSolid,Err)
    604. 

  604.   CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err)
    605. 

  605.   CALL CMISSEquationsSet_CreateStart(EquationSetSolidUserNumber,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
    606. 

  606.     & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
    607. 

  607.     & EquationsSetFieldSolidUserNumber,EquationsSetFieldSolid,EquationsSetSolid,Err)
    608. 

  608.   CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err)
    609. 

  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. 

  610.   !Need to store number of compartments, as finite elasticity uses this to calculate the total mass increase for the constiutive law
    611. 

  611.   CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    612. 

  612.      & CMISS_FIELD_VALUES_SET_TYPE,1,0_CMISSIntg,Err)
    613. 

  613.   CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    614. 

  614.      & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
    615. 

  615. 

  616.   !
    617. 

  617.   !================================================================================================================================
    618. 

  618.   !
    619. 

  619.   ! Solid
    620. 

  620. 

  621.   !Create a dependent field with two variables and four components
    622. 

  622.   CALL CMISSField_Initialise(DependentFieldSolid,Err)
    623. 

  623.   !
    624. 

  624.   CALL CMISSField_CreateStart(FieldDependentSolidUserNumber,Region,DependentFieldSolid,Err)
    625. 

  625.   !
    626. 

  626.   CALL CMISSField_TypeSet(DependentFieldSolid,CMISS_FIELD_GENERAL_TYPE,Err)
    627. 

  627.   CALL CMISSField_MeshDecompositionSet(DependentFieldSolid,Decomposition,Err)
    628. 

  628.   CALL CMISSField_GeometricFieldSet(DependentFieldSolid,GeometricFieldSolid,Err)
    629. 

  629.   CALL CMISSField_DependentTypeSet(DependentFieldSolid,CMISS_FIELD_DEPENDENT_TYPE,Err)
    630. 

  630.   !Create 2N+2 number of variables - 2 for solid, 2N for N Darcy compartments
    631. 

  631.   FieldDependentSolidNumberOfVariables=2*Ncompartments+2
    632. 

  632.   CALL CMISSField_NumberOfVariablesSet(DependentFieldSolid,FieldDependentSolidNumberOfVariables,Err)
    633. 

  633.   !create two variables for each compartment
    634. 

  634.   ALLOCATE(VariableTypes(2*Ncompartments+2))
    635. 

  635.   DO num_var=1,Ncompartments+1
    636. 

  636.      VariableTypes(2*num_var-1)=CMISS_FIELD_U_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
    637. 

  637.      VariableTypes(2*num_var)=CMISS_FIELD_DELUDELN_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
    638. 

  638.   ENDDO
    639. 

  639.   CALL CMISSField_VariableTypesSet(DependentFieldSolid,VariableTypes,Err) 
    640. 

  640. !   CALL CMISSField_VariableTypesSet(DependentFieldSolid,(/CMISS_FIELD_U_VARIABLE_TYPE, &
    641. 

  641. !     & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err)
    642. 

  642.     CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    643. 

  643.        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
    644. 

  644.     CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
    645. 

  645.        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
    646. 

  646.   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
    647. 

  647.   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
    648. 

  648.     & FieldDependentSolidNumberOfComponents,Err)
    649. 

  649.   DO icompartment=3,2*Ncompartments+2
    650. 

  650.     CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
    651. 

  651.   ENDDO
    652. 

  652. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    653. 

  653. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    654. 

  654. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    655. 

  655. 

  656.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidDisplMeshComponentNumber,Err)
    657. 

  657.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidDisplMeshComponentNumber,Err)
    658. 

  658.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidDisplMeshComponentNumber,Err)
    659. 

  659.   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4, &
    660. 

  660.     & CMISS_FIELD_NODE_BASED_INTERPOLATION, &
    661. 

  661.     & Err)
    662. 

  662. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
    663. 

  663. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
    664. 

  664.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidLagrMultMeshComponentNumber,Err)
    665. 

  665. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
    666. 

  666. !     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    667. 

  667. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
    668. 

  668. !     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    669. 

  669. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
    670. 

  670. !     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    671. 

  671.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
    672. 

  672.     & SolidDisplMeshComponentNumber, &
    673. 

  673.     & Err)
    674. 

  674.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
    675. 

  675.     & SolidDisplMeshComponentNumber, &
    676. 

  676.     & Err)
    677. 

  677.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
    678. 

  678.     & SolidDisplMeshComponentNumber, &
    679. 

  679.     & Err)
    680. 

  680.   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
    681. 

  681.     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    682. 

  682. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
    683. 

  683. !     & CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
    684. 

  684. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
    685. 

  685.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
    686. 

  686.     & SolidLagrMultMeshComponentNumber, &
    687. 

  687.     & Err)
    688. 

  688.   !loop over the number of compartments
    689. 

  689.   DO icompartment=3,2*Ncompartments+2
    690. 

  690. !     CALL CMISSField_DimensionSet(DependentFieldSolid,VariableTypes(icompartment), &
    691. 

  691. !        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
    692. 

  692.     !CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
    693. 

  693.     DO componentnum=1,FieldDependentFluidNumberOfComponents-1
    694. 

  694.     !set dimension type
    695. 

  695. !     CALL CMISSField_DimensionSet(DependentField,VariableTypes(icompartment), &
    696. 

  696. !        & CMISS_FIELD_SCALAR_DIMENSION_TYPE,Err)
    697. 

  697.       CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, &
    698. 

  698.        & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    699. 

  699.       CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, & 
    700. 

  700.          & DarcyVelMeshComponentNumber,Err)
    701. 

  701.     ENDDO
    702. 

  702.       CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment), &
    703. 

  703.        & FieldDependentFluidNumberOfComponents, &
    704. 

  704.        & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    705. 

  705. !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
    706. 

  706. !       & FieldDependentFluidNumberOfComponents,MESH_COMPONENT_NUMBER_PRESSURE,Err)
    707. 

  707.     CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
    708. 

  708.       & FieldDependentFluidNumberOfComponents,DarcyMassIncreaseMeshComponentNumber,Err)
    709. 

  709.     
    710. 

  710.   ENDDO
    711. 

  711. 

  712. !   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
    713. 

  713. !   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
    714. 

  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. 

  715. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    716. 

  716. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    717. 

  717. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    718. 

  718. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
    719. 

  719. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    720. 

  720. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    721. 

  721. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    722. 

  722. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
    723. 

  723. 

  724.   !
    725. 

  725.   CALL CMISSField_ScalingTypeSet(DependentFieldSolid,CMISS_FIELD_UNIT_SCALING,Err)
    726. 

  726. 

  727.   CALL CMISSField_CreateFinish(DependentFieldSolid,Err)
    728. 

  728. 

  729.   !
    730. 

  730.   !================================================================================================================================
    731. 

  731.   !
    732. 

  732. 

  733.   CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentSolidUserNumber,DependentFieldSolid,Err)
    734. 

  734.   CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err)
    735. 

  735. 

  736.   CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialUserNumber,MaterialFieldSolid,Err)  
    737. 

  737.   CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err)
    738. 

  738. 

  739.   !
    740. 

  740.   !================================================================================================================================
    741. 

  741.   !
    742. 

  742.   DO icompartment = 1,Ncompartments
    743. 

  743.     CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy(icompartment),FieldDependentSolidUserNumber,&
    744. 

  744.       & DependentFieldSolid,Err)
    745. 

  745.     CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy(icompartment),Err)
    746. 

  746.   ENDDO
    747. 

  747. 

  748.   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS+1
    749. 

  749.     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    750. 

  750.       & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
    751. 

  751.     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    752. 

  752.       & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
    753. 

  753.     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U2_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    754. 

  754.       & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
    755. 

  755.     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U3_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    756. 

  756.       & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
    757. 

  757.   ENDDO
    758. 

  758.   !
    759. 

  759.   !================================================================================================================================
    760. 

  760.   !
    761. 

  761.   ALLOCATE(CouplingCoeffs(Ncompartments,Ncompartments))
    762. 

  762.   IF(Ncompartments==2)THEN
    763. 

  763.     CouplingCoeffs(1,1)=0.0E-01_CMISSDP
    764. 

  764. !     CouplingCoeffs(1,2)=-1.0E-04_CMISSDP
    765. 

  765. !     CouplingCoeffs(2,1)=-1.0E-04_CMISSDP
    766. 

  766.     CouplingCoeffs(1,2)=0.0E-01_CMISSDP
    767. 

  767.     CouplingCoeffs(2,1)=0.0E-01_CMISSDP
    768. 

  768.     CouplingCoeffs(2,2)=0.0E-01_CMISSDP
    769. 

  769.   ELSE IF(Ncompartments==3)THEN
    770. 

  770.     CouplingCoeffs(1,1)=1.0E-02_CMISSDP
    771. 

  771.     CouplingCoeffs(1,2)=1.0E-02_CMISSDP
    772. 

  772.     CouplingCoeffs(1,3)=0.0E-02_CMISSDP
    773. 

  773.     CouplingCoeffs(2,1)=1.0E-02_CMISSDP
    774. 

  774.     CouplingCoeffs(2,2)=2.0E-02_CMISSDP
    775. 

  775.     CouplingCoeffs(2,3)=1.0E-02_CMISSDP
    776. 

  776.     CouplingCoeffs(3,1)=0.0E-02_CMISSDP
    777. 

  777.     CouplingCoeffs(3,2)=1.0E-02_CMISSDP
    778. 

  778.     CouplingCoeffs(3,3)=1.0E-02_CMISSDP
    779. 

  779.   ELSE IF(Ncompartments==4)THEN
    780. 

  780.     CouplingCoeffs(1,1)=0.0E-02_CMISSDP
    781. 

  781.     CouplingCoeffs(1,2)=0.0E-02_CMISSDP
    782. 

  782.     CouplingCoeffs(1,3)=0.0E-02_CMISSDP
    783. 

  783.     CouplingCoeffs(1,4)=0.0E-02_CMISSDP
    784. 

  784.     CouplingCoeffs(2,1)=0.0E-02_CMISSDP
    785. 

  785.     CouplingCoeffs(2,2)=0.0E-02_CMISSDP
    786. 

  786.     CouplingCoeffs(2,3)=0.0E-02_CMISSDP
    787. 

  787.     CouplingCoeffs(2,4)=0.0E-02_CMISSDP
    788. 

  788.     CouplingCoeffs(3,1)=0.0E-02_CMISSDP
    789. 

  789.     CouplingCoeffs(3,2)=0.0E-02_CMISSDP
    790. 

  790.     CouplingCoeffs(3,3)=0.0E-02_CMISSDP
    791. 

  791.     CouplingCoeffs(3,4)=0.0E-02_CMISSDP
    792. 

  792.     CouplingCoeffs(4,1)=0.0E-02_CMISSDP
    793. 

  793.     CouplingCoeffs(4,2)=0.0E-02_CMISSDP
    794. 

  794.     CouplingCoeffs(4,3)=0.0E-02_CMISSDP
    795. 

  795.     CouplingCoeffs(4,4)=0.0E-02_CMISSDP
    796. 

  796.   ELSE IF(Ncompartments==5)THEN
    797. 

  797.     CouplingCoeffs(1,1)=0.0E-02_CMISSDP
    798. 

  798.     CouplingCoeffs(1,2)=0.0E-02_CMISSDP
    799. 

  799.     CouplingCoeffs(1,3)=0.0E-02_CMISSDP
    800. 

  800.     CouplingCoeffs(1,4)=0.0E-02_CMISSDP
    801. 

  801.     CouplingCoeffs(1,5)=0.0E-02_CMISSDP
    802. 

  802.     CouplingCoeffs(2,1)=0.0E-02_CMISSDP
    803. 

  803.     CouplingCoeffs(2,2)=0.0E-02_CMISSDP
    804. 

  804.     CouplingCoeffs(2,3)=0.0E-02_CMISSDP
    805. 

  805.     CouplingCoeffs(2,4)=0.0E-02_CMISSDP
    806. 

  806.     CouplingCoeffs(2,5)=0.0E-02_CMISSDP
    807. 

  807.     CouplingCoeffs(3,1)=0.0E-02_CMISSDP
    808. 

  808.     CouplingCoeffs(3,2)=0.0E-02_CMISSDP
    809. 

  809.     CouplingCoeffs(3,3)=0.0E-02_CMISSDP
    810. 

  810.     CouplingCoeffs(3,4)=0.0E-02_CMISSDP
    811. 

  811.     CouplingCoeffs(3,5)=0.0E-02_CMISSDP
    812. 

  812.     CouplingCoeffs(4,1)=0.0E-02_CMISSDP
    813. 

  813.     CouplingCoeffs(4,2)=0.0E-02_CMISSDP
    814. 

  814.     CouplingCoeffs(4,3)=0.0E-02_CMISSDP
    815. 

  815.     CouplingCoeffs(4,4)=0.0E-02_CMISSDP
    816. 

  816.     CouplingCoeffs(4,5)=0.0E-02_CMISSDP
    817. 

  817.     CouplingCoeffs(5,1)=0.0E-02_CMISSDP
    818. 

  818.     CouplingCoeffs(5,2)=0.0E-02_CMISSDP
    819. 

  819.     CouplingCoeffs(5,3)=0.0E-02_CMISSDP
    820. 

  820.     CouplingCoeffs(5,4)=0.0E-02_CMISSDP
    821. 

  821.     CouplingCoeffs(5,5)=0.0E-02_CMISSDP
    822. 

  822.   ELSE
    823. 

  823.     write(*,*) "Can't initialise coupling coefficients array."
    824. 

  824.   ENDIF
    825. 

  825.   !Define the material parameters for each compartments' constitutive law (for determining pressure)
    826. 

  826.   Nparams=3
    827. 

  827.   ALLOCATE(ConstitutiveParams(Ncompartments,Nparams))
    828. 

  828.   IF(Ncompartments==2)THEN
    829. 

  829.     ConstitutiveParams(1,1)=10.0E-01_CMISSDP
    830. 

  830.     ConstitutiveParams(1,2)=10.0E-01_CMISSDP
    831. 

  831.     ConstitutiveParams(1,3)=10.0E-01_CMISSDP
    832. 

  832.     ConstitutiveParams(2,1)=10.0E-01_CMISSDP
    833. 

  833.     ConstitutiveParams(2,2)=10.0E-01_CMISSDP
    834. 

  834.     ConstitutiveParams(2,3)=10.0E-01_CMISSDP
    835. 

  835.   ELSE IF(Ncompartments==3)THEN
    836. 

  836.     ConstitutiveParams(1,1)=1.0E-02_CMISSDP
    837. 

  837.     ConstitutiveParams(1,2)=1.0E-02_CMISSDP
    838. 

  838.     ConstitutiveParams(1,3)=0.0E-02_CMISSDP
    839. 

  839.     ConstitutiveParams(2,1)=1.0E-02_CMISSDP
    840. 

  840.     ConstitutiveParams(2,2)=2.0E-02_CMISSDP
    841. 

  841.     ConstitutiveParams(2,3)=1.0E-02_CMISSDP
    842. 

  842.     ConstitutiveParams(3,1)=0.0E-02_CMISSDP
    843. 

  843.     ConstitutiveParams(3,2)=1.0E-02_CMISSDP
    844. 

  844.     ConstitutiveParams(3,3)=1.0E-02_CMISSDP
    845. 

  845.   ELSE IF(Ncompartments==4)THEN
    846. 

  846.     ConstitutiveParams(1,1)=0.0E-02_CMISSDP
    847. 

  847.     ConstitutiveParams(1,2)=0.0E-02_CMISSDP
    848. 

  848.     ConstitutiveParams(1,3)=0.0E-02_CMISSDP
    849. 

  849.     ConstitutiveParams(2,1)=0.0E-02_CMISSDP
    850. 

  850.     ConstitutiveParams(2,2)=0.0E-02_CMISSDP
    851. 

  851.     ConstitutiveParams(2,3)=0.0E-02_CMISSDP
    852. 

  852.     ConstitutiveParams(3,1)=0.0E-02_CMISSDP
    853. 

  853.     ConstitutiveParams(3,2)=0.0E-02_CMISSDP
    854. 

  854.     ConstitutiveParams(3,3)=0.0E-02_CMISSDP
    855. 

  855.     ConstitutiveParams(4,1)=0.0E-02_CMISSDP
    856. 

  856.     ConstitutiveParams(4,2)=0.0E-02_CMISSDP
    857. 

  857.     ConstitutiveParams(4,3)=0.0E-02_CMISSDP
    858. 

  858.   ELSE IF(Ncompartments==5)THEN
    859. 

  859.     ConstitutiveParams(1,1)=0.0E-02_CMISSDP
    860. 

  860.     ConstitutiveParams(1,2)=0.0E-02_CMISSDP
    861. 

  861.     ConstitutiveParams(1,3)=0.0E-02_CMISSDP
    862. 

  862.     ConstitutiveParams(2,1)=0.0E-02_CMISSDP
    863. 

  863.     ConstitutiveParams(2,2)=0.0E-02_CMISSDP
    864. 

  864.     ConstitutiveParams(2,3)=0.0E-02_CMISSDP
    865. 

  865.     ConstitutiveParams(3,1)=0.0E-02_CMISSDP
    866. 

  866.     ConstitutiveParams(3,2)=0.0E-02_CMISSDP
    867. 

  867.     ConstitutiveParams(3,3)=0.0E-02_CMISSDP
    868. 

  868.     ConstitutiveParams(4,1)=0.0E-02_CMISSDP
    869. 

  869.     ConstitutiveParams(4,2)=0.0E-02_CMISSDP
    870. 

  870.     ConstitutiveParams(4,3)=0.0E-02_CMISSDP
    871. 

  871.     ConstitutiveParams(5,1)=0.0E-02_CMISSDP
    872. 

  872.     ConstitutiveParams(5,2)=0.0E-02_CMISSDP
    873. 

  873.     ConstitutiveParams(5,3)=0.0E-02_CMISSDP
    874. 

  874.   ELSE
    875. 

  875.     write(*,*) "Can't initialise constitutive parameters array."
    876. 

  876.   ENDIF
    877. 

  877.   !MATERIALS FIELDS
    878. 

  878.   !Auto-created field contains a U variable type to store the diffusion coefficient(s)
    879. 

  879.   !It also contains a V variable type to store the coupling coefficients 
    880. 

  880.   DO icompartment = 1,Ncompartments
    881. 

  881.     MaterialsFieldUserNumberDarcy = 400+icompartment
    882. 

  882.     CALL CMISSField_Initialise(MaterialsFieldDarcy(icompartment),Err)
    883. 

  883.     CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy(icompartment),MaterialsFieldUserNumberDarcy,&
    884. 

  884.          & MaterialsFieldDarcy(icompartment),Err)
    885. 

  885.     CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy(icompartment),Err)
    886. 

  886.   END DO
    887. 

  887.   DO icompartment = 1,Ncompartments
    888. 

  888.     CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    889. 

  889.       & CMISS_FIELD_VALUES_SET_TYPE, &
    890. 

  890.       & MaterialsFieldUserNumberDarcyPorosity,POROSITY_PARAM_DARCY,Err)
    891. 

  891.     CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    892. 

  892.       & CMISS_FIELD_VALUES_SET_TYPE, &
    893. 

  893.       & MaterialsFieldUserNumberDarcyPermOverVis,PERM_OVER_VIS_PARAM_DARCY,Err)
    894. 

  894.   END DO
    895. 

  895.   DO icompartment = 1, Ncompartments
    896. 

  896.     DO COMPONENT_NUMBER=1, Ncompartments
    897. 

  897.       CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
    898. 

  898.          & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
    899. 

  899. !         CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
    900. 

  900. !           & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
    901. 

  901.     END DO
    902. 

  902.   END DO
    903. 

  903.   DO icompartment = 1, Ncompartments
    904. 

  904.     DO COMPONENT_NUMBER=1,Nparams
    905. 

  905.       CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U1_VARIABLE_TYPE, &
    906. 

  906.          & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,ConstitutiveParams(icompartment,COMPONENT_NUMBER),Err)
    907. 

  907. !         CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
    908. 

  908. !           & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
    909. 

  909.     END DO
    910. 

  910.   END DO
    911. 

  911. 

  912.   !
    913. 

  913.   !================================================================================================================================
    914. 

  914.   !
    915. 

  915. 

  916.   !EQUATIONS SET EQUATIONS
    917. 

  917. 

  918.   !Darcy
    919. 

  919.   DO icompartment=1,Ncompartments
    920. 

  920.     !Create the equations set equations
    921. 

  921.     CALL CMISSEquations_Initialise(EquationsDarcy(icompartment),Err)
    922. 

  922.     CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy(icompartment),EquationsDarcy(icompartment),Err)
    923. 

  923.     !Set the equations matrices sparsity type
    924. 

  924.     CALL CMISSEquations_SparsityTypeSet(EquationsDarcy(icompartment…
    925.