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

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

Large files are truncated click here to view the full file

    

  1. !> \file
    2. 

  2. !> \author Christian Michler, Adam Reeve, Andrew Cookson
    3. 

  3. !> \brief This is an example program to solve a coupled Finite Elastiticity Darcy equation 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):
    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/IncompressibleElasticityDrivenDarcy/src/IncompressibleElasticityDrivenDarcyExample.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/IncompressibleElasticityDrivenDarcy/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/IncompressibleElasticityDrivenDarcy/build-intel'>Linux GNU Build</a>
    47. 

  47. !!
    48. 

  48. !<
    49. 

  49. 

  50. ! !
    51. 

  51. ! !  This example considers a coupled Finite Elasticity Multi-Compartment Darcy problem
    52. 

  52. ! !
    53. 

  53. 

  54. !> Main program
    55. 

  55. 

  56. PROGRAM FINITEELASTICITYMULTICOMPDARCYEXAMPLE
    57. 

  57. 

  58.   !
    59. 

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

  60.   !
    61. 

  61. 

  62.   !PROGRAM LIBRARIES
    63. 

  63. 

  64.   USE OPENCMISS
    65. 

  65.   USE FLUID_MECHANICS_IO_ROUTINES
    66. 

  66.   USE MPI
    67. 

  67. 

  68. #ifdef WIN32
    69. 

  69.   USE IFQWINCMISS
    70. 

  70. #endif
    71. 

  71. 

  72.   !
    73. 

  73.   !================================================================================================================================
    74. 

  74.   !
    75. 

  75. 

  76.   !PROGRAM VARIABLES AND TYPES
    77. 

  77. 

  78.   IMPLICIT NONE
    79. 

  79. 

  80.   !Test program parameters
    81. 

  81. 

  82.   REAL(CMISSDP), PARAMETER :: Y_DIM=1.0_CMISSDP
    83. 

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

  84.   REAL(CMISSDP), PARAMETER :: Z_DIM=3.0_CMISSDP
    85. 

  85. 

  86.   INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=1
    87. 

  87.   INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=2
    88. 

  88.   INTEGER(CMISSIntg), PARAMETER :: CubicBasisUserNumber=3
    89. 

  89. 

  90.   INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
    91. 

  91.   INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=2
    92. 

  92.   INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=3
    93. 

  93.   INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=4
    94. 

  94.   INTEGER(CMISSIntg), PARAMETER :: GeometricFieldUserNumber=5
    95. 

  95.   INTEGER(CMISSIntg), PARAMETER :: DependentFieldUserNumberDarcy=6
    96. 

  96.   INTEGER(CMISSIntg) :: MaterialsFieldUserNumberDarcy
    97. 

  97.   INTEGER(CMISSIntg) :: EquationsSetUserNumberDarcy
    98. 

  98.   INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=14
    99. 

  99.   INTEGER(CMISSIntg) :: EquationsSetFieldUserNumberDarcy
    100. 

  100.   INTEGER(CMISSIntg) :: icompartment,Ncompartments,num_var,componentnum,Nparams
    101. 

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

  102.   INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2
    103. 

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

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

  105.   INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1
    106. 

  106.   INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPorosity=1
    107. 

  107.   INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPermOverVis=2
    108. 

  108.   INTEGER(CMISSIntg) :: SourceFieldDarcyUserNumber
    109. 

  109. 

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

  111.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
    112. 

  112.   INTEGER(CMISSIntg) :: IndependentFieldDarcyUserNumber
    113. 

  113.   !Program types
    114. 

  114. 

  115.   INTEGER(CMISSIntg) :: NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS,NUMBER_GLOBAL_Z_ELEMENTS
    116. 

  116. 

  117. !   INTEGER(CMISSIntg) :: MPI_IERROR
    118. 

  118.   INTEGER(CMISSIntg) :: NumberOfComputationalNodes,NumberOfDomains,ComputationalNodeNumber
    119. 

  119. 

  120.   TYPE(EXPORT_CONTAINER):: CM
    121. 

  121. 

  122.   !Program variables
    123. 

  123. 

  124.   INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
    125. 

  125.   INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
    126. 

  126.   INTEGER(CMISSIntg) :: RESTART_VALUE
    127. 

  127. 

  128.   INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT
    129. 

  129.   INTEGER(CMISSIntg) :: COMPONENT_NUMBER
    130. 

  130.   INTEGER(CMISSIntg) :: NODE_NUMBER
    131. 

  131.   INTEGER(CMISSIntg) :: ELEMENT_NUMBER
    132. 

  132.   INTEGER(CMISSIntg) :: CONDITION
    133. 

  133. 

  134.   INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY
    135. 

  135.   INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE
    136. 

  136.   INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE
    137. 

  137. 

  138.   REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z
    139. 

  139.   REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2
    140. 

  140.   REAL(CMISSDP) :: GEOMETRY_TOLERANCE
    141. 

  141.   INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID
    142. 

  142.   REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4)
    143. 

  143.   REAL(CMISSDP) :: INITIAL_FIELD_SOLID(4)
    144. 

  144.   REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
    145. 

  145.   REAL(CMISSDP) :: RELATIVE_TOLERANCE
    146. 

  146.   REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
    147. 

  147.   REAL(CMISSDP) :: LINESEARCH_ALPHA
    148. 

  148.   REAL(CMISSDP) :: VALUE
    149. 

  149.   REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY
    150. 

  150. 

  151.   LOGICAL :: EXPORT_FIELD_IO
    152. 

  152.   LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG
    153. 

  153. 

  154.   !CMISS variables
    155. 

  155. 

  156.   !Regions
    157. 

  157.   TYPE(CMISSRegionType) :: Region
    158. 

  158.   TYPE(CMISSRegionType) :: WorldRegion
    159. 

  159.   !Coordinate systems
    160. 

  160.   TYPE(CMISSCoordinateSystemType) :: CoordinateSystem
    161. 

  161.   TYPE(CMISSCoordinateSystemType) :: WorldCoordinateSystem
    162. 

  162.   !Basis
    163. 

  163.   TYPE(CMISSBasisType) :: CubicBasis, QuadraticBasis, LinearBasis, Bases(2)
    164. 

  164.   !Meshes
    165. 

  165.   TYPE(CMISSMeshType) :: Mesh
    166. 

  166.   TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
    167. 

  167.   !Decompositions
    168. 

  168.   TYPE(CMISSDecompositionType) :: Decomposition
    169. 

  169.   !Fields
    170. 

  170.   TYPE(CMISSFieldsType) :: Fields
    171. 

  171.   !Field types
    172. 

  172.   TYPE(CMISSFieldType) :: GeometricField
    173. 

  173.   TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: MaterialsFieldDarcy
    174. 

  174.   TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: EquationsSetFieldDarcy
    175. 

  175.   TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: SourceFieldDarcy
    176. 

  176.   TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: IndependentFieldDarcy
    177. 

  177.   !Boundary conditions
    178. 

  178.   TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy
    179. 

  179.   !Equations sets
    180. 

  180.   TYPE(CMISSEquationsSetType), ALLOCATABLE, DIMENSION(:) :: EquationsSetDarcy
    181. 

  181.   !Equations
    182. 

  182.   TYPE(CMISSEquationsType), ALLOCATABLE, DIMENSION(:) :: EquationsDarcy
    183. 

  183.   !Problems
    184. 

  184.   TYPE(CMISSProblemType) :: Problem
    185. 

  185.   !Control loops
    186. 

  186.   TYPE(CMISSControlLoopType) :: ControlLoop
    187. 

  187.   !Solvers
    188. 

  188.   TYPE(CMISSSolverType) :: DynamicSolverDarcy
    189. 

  189.   TYPE(CMISSSolverType) :: LinearSolverDarcy
    190. 

  190. !   TYPE(CMISSSolverType) :: LinearSolverSolid
    191. 

  191.   !Solver equations
    192. 

  192.   TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy
    193. 

  193.   !Other variables
    194. 

  194.   INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face1Nodes(:),Face2Nodes(:)
    195. 

  195.   INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face3Nodes(:),Face4Nodes(:)
    196. 

  196.   INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face5Nodes(:),Face6Nodes(:)
    197. 

  197.   INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face7Nodes(:),Face8Nodes(:)
    198. 

  198.   INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face9Nodes(:),Face10Nodes(:)
    199. 

  199.   INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face11Nodes(:),Face12Nodes(:)
    200. 

  200.   INTEGER(CMISSIntg) :: FaceXi(6)
    201. 

  201.   INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
    202. 

  202.   REAL(CMISSDP) :: XCoord,YCoord,ZCoord
    203. 

  203.   LOGICAL :: X_FIXED,Y_FIXED !,X_OKAY,Y_OKAY
    204. 

  204. 

  205. #ifdef WIN32
    206. 

  206.   !Quickwin type
    207. 

  207.   LOGICAL :: QUICKWIN_STATUS=.FALSE.
    208. 

  208.   TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
    209. 

  209. #endif
    210. 

  210. 

  211.   !Generic CMISS variables
    212. 

  212. 

  213.   INTEGER(CMISSIntg) :: EquationsSetIndex
    214. 

  214.   INTEGER(CMISSIntg) :: Err
    215. 

  215.   !Array containing the field variable types that will be used (for ease of incorporating inside a loop)
    216. 

  216.   INTEGER(CMISSIntg), ALLOCATABLE, DIMENSION(:) :: VariableTypes
    217. 

  217.   REAL(CMISSDP), ALLOCATABLE, DIMENSION(:,:) :: CouplingCoeffs,ConstitutiveParams
    218. 

  218. 

  219.   INTEGER(CMISSIntg) :: DIAG_LEVEL_LIST(5)
    220. 

  220. !   CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(8) !,TIMING_ROUTINE_LIST(1)
    221. 

  221.   CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(1) !,TIMING_ROUTINE_LIST(1)
    222. 

  222. 

  223.   !
    224. 

  224.   !--------------------------------------------------------------------------------------------------------------------------------
    225. 

  225.   !
    226. 

  226. 

  227.   !Program variables and types (finite elasticity part)
    228. 

  228. 

  229.   !Test program parameters
    230. 

  230. 

  231.   INTEGER(CMISSIntg) :: SolidMeshComponenetNumber
    232. 

  232. 

  233.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidUserNumber=51
    234. 

  234.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfVariables=1
    235. 

  235.   INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfComponents=3
    236. 

  236. 

  237.   INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidUserNumber=52
    238. 

  238.   INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfVariables=1
    239. 

  239.   INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfComponents=3
    240. 

  240. 

  241.   INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidUserNumber=53
    242. 

  242.   INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfVariables=1
    243. 

  243.   INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfComponents=3
    244. 

  244. 

  245.   INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidUserNumber=54
    246. 

  246.   INTEGER(CMISSIntg) :: FieldDependentSolidNumberOfVariables
    247. 

  247.   INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4
    248. 

  248.   INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4  !(u,v,w,m)
    249. 

  249. 

  250.   INTEGER(CMISSIntg), PARAMETER :: EquationSetSolidUserNumber=55
    251. 

  251.   INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldSolidUserNumber=25
    252. 

  252. 

  253.   INTEGER(CMISSIntg), PARAMETER :: SolidDisplMeshComponentNumber=1
    254. 

  254.   INTEGER(CMISSIntg), PARAMETER :: SolidLagrMultMeshComponentNumber=2
    255. 

  255.   INTEGER(CMISSIntg), PARAMETER :: SolidGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
    256. 

  256. 

  257.   INTEGER(CMISSIntg), PARAMETER :: DarcyVelMeshComponentNumber=SolidLagrMultMeshComponentNumber
    258. 

  258.   INTEGER(CMISSIntg), PARAMETER :: DarcyMassIncreaseMeshComponentNumber=SolidLagrMultMeshComponentNumber
    259. 

  259. !   INTEGER(CMISSIntg), PARAMETER :: DarcyGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
    260. 

  260. 

  261.   INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=32
    262. 

  262.   !Program types
    263. 

  263.   !Program variables
    264. 

  264. 

  265.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME
    266. 

  266.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME
    267. 

  267.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA
    268. 

  268.   REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
    269. 

  269. 

  270.   !CMISS variables
    271. 

  271. 

  272. !   TYPE(CMISSBasisType) :: BasisSolid
    273. 

  273.   TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsSolid
    274. 

  274.   TYPE(CMISSEquationsType) :: EquationsSolid
    275. 

  275.   TYPE(CMISSEquationsSetType) :: EquationsSetSolid
    276. 

  276.   TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid,DependentFieldSolid,EquationsSetFieldSolid
    277. 

  277.   TYPE(CMISSSolverType) :: SolverSolid
    278. 

  278.   TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid
    279. 

  279.   TYPE(CMISSMeshElementsType) :: MeshElementsSolid
    280. 

  280. 

  281.   !End - Program variables and types (finite elasticity part)
    282. 

  282. 

  283.   INTEGER(CMISSIntg) :: dummy_counter
    284. 

  284. 

  285.   !
    286. 

  286.   !--------------------------------------------------------------------------------------------------------------------------------
    287. 

  287.   !
    288. 

  288. 

  289. 

  290. #ifdef WIN32
    291. 

  291.   !Initialise QuickWin
    292. 

  292.   QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
    293. 

  293.   QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
    294. 

  294.   QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
    295. 

  295.   !Set the window parameters
    296. 

  296.   QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
    297. 

  297.   !If attempt fails set with system estimated values
    298. 

  298.   IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
    299. 

  299. #endif
    300. 

  300. 

  301.   !
    302. 

  302.   !================================================================================================================================
    303. 

  303.   !
    304. 

  304. 

  305.   !INITIALISE OPENCMISS
    306. 

  306. 

  307.   CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
    308. 

  308. 

  309.   !CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
    310. 

  310. 

  311.   !
    312. 

  312.   !================================================================================================================================
    313. 

  313.   !
    314. 

  314. 

  315.   !PROBLEM CONTROL PANEL
    316. 

  316.   NUMBER_GLOBAL_X_ELEMENTS=1
    317. 

  317.   NUMBER_GLOBAL_Y_ELEMENTS=1
    318. 

  318.   NUMBER_GLOBAL_Z_ELEMENTS=3
    319. 

  319. 

  320.   IF(NUMBER_GLOBAL_Z_ELEMENTS==0)THEN
    321. 

  321.     NUMBER_OF_DIMENSIONS=2
    322. 

  322.   ELSE
    323. 

  323.     NUMBER_OF_DIMENSIONS=3
    324. 

  324.   ENDIF
    325. 

  325.   !PROBLEM CONTROL PANEL
    326. 

  326. 

  327. !   BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION
    328. 

  328.   BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION
    329. 

  329.   !Set geometric tolerance
    330. 

  330.   GEOMETRY_TOLERANCE = 1.0E-12_CMISSDP
    331. 

  331.   !Set initial values
    332. 

  332.   INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
    333. 

  333.   INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
    334. 

  334.   INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
    335. 

  335.   INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
    336. 

  336.   !Set material parameters
    337. 

  337.   POROSITY_PARAM_DARCY=0.1_CMISSDP
    338. 

  338.   PERM_OVER_VIS_PARAM_DARCY=1.0_CMISSDP
    339. 

  339.   !Set output parameter
    340. 

  340.   !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
    341. 

  341.   DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_MATRIX_OUTPUT
    342. 

  342.   LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_NO_OUTPUT
    343. 

  343.   !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
    344. 

  344.   EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
    345. 

  345. 

  346.   !Set time parameter
    347. 

  347.   DYNAMIC_SOLVER_DARCY_START_TIME=1.0E-3_CMISSDP
    348. 

  348.   DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP
    349. 

  349.   DYNAMIC_SOLVER_DARCY_STOP_TIME=2_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
    350. 

  350.   DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
    351. 

  351.   !Set result output parameter
    352. 

  352.   DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
    353. 

  353.   !Set solver parameters
    354. 

  354.   LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
    355. 

  355.   RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
    356. 

  356.   ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
    357. 

  357.   DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
    358. 

  358.   MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
    359. 

  359.   RESTART_VALUE=30_CMISSIntg !default: 30
    360. 

  360.   LINESEARCH_ALPHA=1.0_CMISSDP
    361. 

  361. 

  362. !   !Import cmHeart mesh information
    363. 

  363. !   CALL FLUID_MECHANICS_IO_READ_CMHEART(CM,Err)
    364. 

  364. !   BASIS_NUMBER_GEOMETRY=CM%ID_M
    365. 

  365. !   BASIS_NUMBER_VELOCITY=CM%ID_V
    366. 

  366. !   BASIS_NUMBER_PRESSURE=CM%ID_P
    367. 

  367. !   NUMBER_OF_DIMENSIONS=CM%D
    368. 

  368. !   BASIS_TYPE=CM%IT_T
    369. 

  369. !   BASIS_XI_INTERPOLATION_GEOMETRY=CM%IT_M
    370. 

  370. !   BASIS_XI_INTERPOLATION_VELOCITY=CM%IT_V
    371. 

  371. !   BASIS_XI_INTERPOLATION_PRESSURE=CM%IT_P
    372. 

  372. !   BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION
    373. 

  373. !   NUMBER_OF_NODES_GEOMETRY=CM%N_M
    374. 

  374. !   NUMBER_OF_NODES_VELOCITY=CM%N_V
    375. 

  375. !   NUMBER_OF_NODES_PRESSURE=CM%N_P
    376. 

  376. !   TOTAL_NUMBER_OF_NODES=CM%N_T
    377. 

  377. !   TOTAL_NUMBER_OF_ELEMENTS=CM%E_T
    378. 

  378. !   NUMBER_OF_ELEMENT_NODES_GEOMETRY=CM%EN_M
    379. 

  379. !   NUMBER_OF_ELEMENT_NODES_VELOCITY=CM%EN_V
    380. 

  380. !   NUMBER_OF_ELEMENT_NODES_PRESSURE=CM%EN_P
    381. 

  381. !   !Set domain dimensions
    382. 

  382. !   DOMAIN_X1 = -5.0_CMISSDP
    383. 

  383. !   DOMAIN_X2 =  5.0_CMISSDP
    384. 

  384. !   DOMAIN_Y1 = -5.0_CMISSDP
    385. 

  385. !   DOMAIN_Y2 =  5.0_CMISSDP
    386. 

  386. !   DOMAIN_Z1 = -5.0_CMISSDP
    387. 

  387. !   DOMAIN_Z2 =  5.0_CMISSDP
    388. 

  388. !   !Set geometric tolerance
    389. 

  389. !   GEOMETRY_TOLERANCE = 1.0E-12_CMISSDP
    390. 

  390. !   !Set initial values
    391. 

  391. !   INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
    392. 

  392. !   INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
    393. 

  393. !   INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
    394. 

  394. !   INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
    395. 

  395. !   INITIAL_FIELD_MAT_PROPERTIES(1)=0.0_CMISSDP
    396. 

  396. !   INITIAL_FIELD_MAT_PROPERTIES(2)=0.0_CMISSDP
    397. 

  397. !   INITIAL_FIELD_MAT_PROPERTIES(3)=0.0_CMISSDP
    398. 

  398. ! !   INITIAL_FIELD_SOLID(1)=1.0_CMISSDP
    399. 

  399. ! !   INITIAL_FIELD_SOLID(2)=1.0_CMISSDP
    400. 

  400. ! !   INITIAL_FIELD_SOLID(3)=1.0_CMISSDP
    401. 

  401. ! !   INITIAL_FIELD_SOLID(4)=1.0_CMISSDP
    402. 

  402. !   !Set material parameters
    403. 

  403. !   POROSITY_PARAM_DARCY=0.1_CMISSDP
    404. 

  404. !   PERM_OVER_VIS_PARAM_DARCY=1.0e-1_CMISSDP
    405. 

  405. !   POROSITY_PARAM_MAT_PROPERTIES=POROSITY_PARAM_DARCY
    406. 

  406. !   PERM_OVER_VIS_PARAM_MAT_PROPERTIES=PERM_OVER_VIS_PARAM_DARCY
    407. 

  407. !   !Set number of Gauss points (Mind that also material field may be interpolated)
    408. 

  408. !   BASIS_XI_GAUSS_GEOMETRY=3 !4
    409. 

  409. !   BASIS_XI_GAUSS_VELOCITY=3 !4
    410. 

  410. !   BASIS_XI_GAUSS_PRESSURE=3 !4
    411. 

  411. !   !Set output parameter
    412. 

  412. !   !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
    413. 

  413. !   LINEAR_SOLVER_MAT_PROPERTIES_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
    414. 

  414. !   DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
    415. 

  415. !   LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_SOLVER_OUTPUT
    416. 

  416. !   !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
    417. 

  417. !   EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
    418. 

  418. !   EQUATIONS_MAT_PROPERTIES_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
    419. 

  419. ! 
    420. 

  420. !   !Set time parameter
    421. 

  421. !   DYNAMIC_SOLVER_DARCY_START_TIME=0.0_CMISSDP
    422. 

  422. ! !   DYNAMIC_SOLVER_DARCY_STOP_TIME=0.03_CMISSDP
    423. 

  423. !   DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-2_CMISSDP
    424. 

  424. !   DYNAMIC_SOLVER_DARCY_STOP_TIME=1_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
    425. 

  425. !   DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
    426. 

  426. !   !Set result output parameter
    427. 

  427. !   DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
    428. 

  428. !   !Set solver parameters
    429. 

  429. !   LINEAR_SOLVER_MAT_PROPERTIES_DIRECT_FLAG=.TRUE.
    430. 

  430. !   LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
    431. 

  431. !   RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
    432. 

  432. !   ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
    433. 

  433. !   DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
    434. 

  434. !   MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
    435. 

  435. !   RESTART_VALUE=30_CMISSIntg !default: 30
    436. 

  436. !   LINESEARCH_ALPHA=1.0_CMISSDP
    437. 

  437. 

  438.   icompartment =1_CMISSIntg
    439. 

  439.   Ncompartments=2_CMISSIntg
    440. 

  440.   !
    441. 

  441.   !================================================================================================================================
    442. 

  442.   !
    443. 

  443. 

  444. ! !   !Set diagnostics
    445. 

  445. ! ! 
    446. 

  446. ! !   DIAG_LEVEL_LIST(1)=1
    447. 

  447. ! !   DIAG_LEVEL_LIST(2)=2
    448. 

  448. ! !   DIAG_LEVEL_LIST(3)=3
    449. 

  449. ! !   DIAG_LEVEL_LIST(4)=4
    450. 

  450. ! !   DIAG_LEVEL_LIST(5)=5
    451. 

  451. ! ! 
    452. 

  452. ! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_FINITE_ELEMENT_CALCULATE"
    453. 

  453. ! ! !   DIAG_ROUTINE_LIST(2)="DARCY_EQUATION_PRE_SOLVE_STORE_REFERENCE_DATA"
    454. 

  454. ! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_GET_SOLID_DISPLACEMENT"
    455. 

  455. ! ! !   DIAG_ROUTINE_LIST(2)="DARCY_EQUATION_PRE_SOLVE_ALE_UPDATE_MESH"
    456. 

  456. ! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_GET_SOLID_DISPLACEMENT"
    457. 

  457. ! ! !   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_UPDATE_BOUNDARY_CONDITIONS"
    458. 

  458. ! !   DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
    459. 

  459. ! ! !   DIAG_ROUTINE_LIST(2)="FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR"
    460. 

  460. ! ! !   DIAG_ROUTINE_LIST(3)="EVALUATE_CHAPELLE_PIOLA_TENSOR_ADDITION"
    461. 

  461. ! ! !   DIAG_ROUTINE_LIST(5)="DARCY_EQUATION_PRE_SOLVE_MAT_PROPERTIES"
    462. 

  462. ! ! !   DIAG_ROUTINE_LIST(6)="DATA_FITTING_FINITE_ELEMENT_CALCULATE"
    463. 

  463. ! ! !   DIAG_ROUTINE_LIST(7)="FINITE_ELASTICITY_FINITE_ELEMENT_JACOBIAN_EVALUATE"
    464. 

  464. ! ! !   DIAG_ROUTINE_LIST(8)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
    465. 

  465. ! ! !   DIAG_ROUTINE_LIST(1)="PROBLEM_SOLVER_EQUATIONS_SOLVE"
    466. 

  466. ! ! !   DIAG_ROUTINE_LIST(1)="SOLVER_NEWTON_SOLVE"
    467. 

  467. ! ! !   DIAG_ROUTINE_LIST(2)="SOLVER_NEWTON_LINESEARCH_SOLVE"
    468. 

  468. ! ! !   DIAG_ROUTINE_LIST(1)="SOLVER_SOLUTION_UPDATE"
    469. 

  469. ! ! !   DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
    470. 

  470. 

  471. 

  472. !Set diagnostics
    473. 

  473. 

  474.   DIAG_LEVEL_LIST(1)=1
    475. 

  475.   DIAG_LEVEL_LIST(2)=2
    476. 

  476.   DIAG_LEVEL_LIST(3)=3
    477. 

  477.   DIAG_LEVEL_LIST(4)=4
    478. 

  478.   DIAG_LEVEL_LIST(5)=5
    479. 

  479. 

  480.   !DIAG_ROUTINE_LIST(1)="WRITE_IP_INFO"
    481. 

  481. !   DIAG_ROUTINE_LIST(2)="FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR"
    482. 

  482.   DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
    483. 

  483. 

  484.   !CMISS_ALL_DIAG_TYPE/CMISS_IN_DIAG_TYPE/CMISS_FROM_DIAG_TYPE
    485. 

  485.   CALL CMISSDiagnosticsSetOn(CMISS_IN_DIAG_TYPE,DIAG_LEVEL_LIST,"Diagnostics",DIAG_ROUTINE_LIST,Err)
    486. 

  486. 

  487.   !CMISS_ALL_TIMING_TYPE/CMISS_IN_TIMING_TYPE/CMISS_FROM_TIMING_TYPE
    488. 

  488.   !TIMING_ROUTINE_LIST(1)="PROBLEM_FINITE_ELEMENT_CALCULATE"
    489. 

  489.   !CALL TIMING_SET_ON(IN_TIMING_TYPE,.TRUE.,"",TIMING_ROUTINE_LIST,ERR,ERROR,*999)
    490. 

  490. 

  491. !   !CMISS_ALL_DIAG_TYPE/CMISS_IN_DIAG_TYPE/CMISS_FROM_DIAG_TYPE
    492. 

  492. !    CALL CMISSDiagnosticsSetOn(CMISS_IN_DIAG_TYPE,DIAG_LEVEL_LIST,"Diagnostics",DIAG_ROUTINE_LIST,Err)
    493. 

  493. 

  494.   !CMISS_ALL_TIMING_TYPE/CMISS_IN_TIMING_TYPE/CMISS_FROM_TIMING_TYPE
    495. 

  495.   !TIMING_ROUTINE_LIST(1)="PROBLEM_FINITE_ELEMENT_CALCULATE"
    496. 

  496.   !CALL TIMING_SET_ON(IN_TIMING_TYPE,.TRUE.,"",TIMING_ROUTINE_LIST,ERR,ERROR,*999)
    497. 

  497.   ALLOCATE (EquationsSetDarcy(Ncompartments))
    498. 

  498.   ALLOCATE (EquationsSetFieldDarcy(Ncompartments))
    499. 

  499.   ALLOCATE (MaterialsFieldDarcy(Ncompartments))
    500. 

  500.   ALLOCATE (EquationsDarcy(Ncompartments))
    501. 

  501.   ALLOCATE (SourceFieldDarcy(Ncompartments))
    502. 

  502.   ALLOCATE (IndependentFieldDarcy(Ncompartments))
    503. 

  503.   !
    504. 

  504.   !================================================================================================================================
    505. 

  505.   !
    506. 

  506. 

  507.   !Get the number of computational nodes and this computational node number
    508. 

  508.   CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
    509. 

  509.   CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
    510. 

  510. 

  511.   NumberOfDomains = NumberOfComputationalNodes
    512. 

  512.   write(*,*) "NumberOfDomains = ",NumberOfDomains
    513. 

  513. 

  514. 

  515.   !
    516. 

  516.   !================================================================================================================================
    517. 

  517.   !
    518. 

  518. 

  519.   !COORDINATE SYSTEM
    520. 

  520. 

  521.   CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
    522. 

  522.   CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
    523. 

  523.   CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NUMBER_OF_DIMENSIONS,Err)
    524. 

  524.   CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
    525. 

  525. 

  526.   !
    527. 

  527.   !================================================================================================================================
    528. 

  528.   !
    529. 

  529. 

  530.   !REGION
    531. 

  531.   !For a volume-coupled problem, solid and fluid are based in the same region
    532. 

  532. 

  533.   CALL CMISSRegion_Initialise(Region,Err)
    534. 

  534.   CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
    535. 

  535.   CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
    536. 

  536.   CALL CMISSRegion_CreateFinish(Region,Err)
    537. 

  537. 

  538.   !
    539. 

  539.   !================================================================================================================================
    540. 

  540.   !
    541. 

  541. 

  542.   !BASES
    543. 

  543.   !Define basis functions
    544. 

  544.   CALL CMISSBasis_Initialise(LinearBasis,Err)
    545. 

  545.   CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
    546. 

  546.   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
    547. 

  547.     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
    548. 

  548.   !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err)
    549. 

  549.   CALL CMISSBasis_CreateFinish(LinearBasis,Err)
    550. 

  550. 

  551.   CALL CMISSBasis_Initialise(QuadraticBasis,Err)
    552. 

  552.   CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
    553. 

  553.   CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
    554. 

  554.     & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
    555. 

  555.   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
    556. 

  556.     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
    557. 

  557.   !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err)
    558. 

  558.   CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
    559. 

  559. 

  560.   CALL CMISSBasis_Initialise(CubicBasis,Err)
    561. 

  561.   CALL CMISSBasis_CreateStart(CubicBasisUserNumber,CubicBasis,Err)
    562. 

  562.   CALL CMISSBasis_InterpolationXiSet(CubicBasis,(/CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION, &
    563. 

  563.     & CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION/),Err)
    564. 

  564.   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(CubicBasis, &
    565. 

  565.     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
    566. 

  566.   !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(CubicBasis,.true.,Err) !Enable 3D interpolation on faces
    567. 

  567.   CALL CMISSBasis_CreateFinish(CubicBasis,Err)
    568. 

  568. 

  569.   !LinearBasis/QuadraticBasis/CubicBasis
    570. 

  570.   Bases(1)=QuadraticBasis
    571. 

  571.   Bases(2)=LinearBasis
    572. 

  572. !   Bases(1)=CubicBasis
    573. 

  573. !   Bases(2)=QuadraticBasis
    574. 

  574. 

  575.   !Start the creation of a generated mesh in the region
    576. 

  576.   CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
    577. 

  577.   CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
    578. 

  578.   CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_REGULAR_MESH_TYPE,Err)
    579. 

  579.   CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,Bases,Err)
    580. 

  580.   IF(NUMBER_GLOBAL_Z_ELEMENTS==0) THEN
    581. 

  581.     CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/X_DIM,Y_DIM/),Err)
    582. 

  582.     CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS/),Err)
    583. 

  583.   ELSE
    584. 

  584.     CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/X_DIM,Y_DIM,Z_DIM/),Err)
    585. 

  585.     CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS, &
    586. 

  586.       & NUMBER_GLOBAL_Z_ELEMENTS/),Err)
    587. 

  587.   ENDIF
    588. 

  588.   CALL CMISSMesh_Initialise(Mesh,Err)
    589. 

  589.   CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
    590. 

  590. 

  591.   !
    592. 

  592.   !================================================================================================================================
    593. 

  593.   !
    594. 

  594. 

  595.   !GEOMETRIC FIELD
    596. 

  596. 

  597.   !Create a decomposition:
    598. 

  598.   !All mesh components (associated with G.Projection / Darcy / solid) share the same decomposition
    599. 

  599.   CALL CMISSDecomposition_Initialise(Decomposition,Err)
    600. 

  600.   CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
    601. 

  601.   !Set the decomposition to be a general decomposition with the specified number of domains
    602. 

  602.   CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
    603. 

  603.   CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
    604. 

  604.   CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.true.,Err)
    605. 

  605.   !Finish the decomposition
    606. 

  606.   CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
    607. 

  607. 

  608.   CALL CMISSField_Initialise(GeometricField,Err)
    609. 

  609.   CALL CMISSField_CreateStart(GeometricFieldUserNumber,Region,GeometricField,Err)
    610. 

  610.   CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
    611. 

  611.   CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
    612. 

  612.   CALL CMISSField_NumberOfVariablesSet(GeometricField,FieldGeometryNumberOfVariables,Err)
    613. 

  613.   CALL CMISSField_NumberOfComponentsSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,Err)  
    614. 

  614.   CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
    615. 

  615.   CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
    616. 

  616.   CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
    617. 

  617.   CALL CMISSField_CreateFinish(GeometricField,Err)
    618. 

  618.   CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricField,Err)
    619. 

  619. 

  620.   !--------------------------------------------------------------------------------------------------------------------------------
    621. 

  621.   ! Solid
    622. 

  622. 

  623.   !Create a decomposition
    624. 

  624. 

  625.   !Create a field to put the geometry (defualt is geometry)
    626. 

  626. 

  627.   SolidMeshComponenetNumber = SolidGeometryMeshComponentNumber
    628. 

  628. 

  629.   CALL CMISSField_Initialise(GeometricFieldSolid,Err)
    630. 

  630.   CALL CMISSField_CreateStart(FieldGeometrySolidUserNumber,Region,GeometricFieldSolid,Err)
    631. 

  631.   CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err)
    632. 

  632.   CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err)
    633. 

  633.   CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometrySolidNumberOfVariables,Err)
    634. 

  634.   CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometrySolidNumberOfComponents,Err)
    635. 

  635.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidMeshComponenetNumber,Err)
    636. 

  636.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidMeshComponenetNumber,Err)
    637. 

  637.   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidMeshComponenetNumber,Err)
    638. 

  638.   CALL CMISSField_CreateFinish(GeometricFieldSolid,Err)
    639. 

  639.   !Set the mesh component to be used by the field components.
    640. 

  640.   CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err)
    641. 

  641. !   !Set the scaling to use
    642. 

  642. !   CALL CMISSField_ScalingTypeSet(GeometricFieldSolid,CMISS_FIELD_NO_SCALING,Err)
    643. 

  643. 

  644.   !Create a fibre field and attach it to the geometric field
    645. 

  645.   CALL CMISSField_Initialise(FibreFieldSolid,Err)
    646. 

  646.   CALL CMISSField_CreateStart(FieldFibreSolidUserNumber,Region,FibreFieldSolid,Err)
    647. 

  647.   CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err)
    648. 

  648.   CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err)
    649. 

  649.   CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err)
    650. 

  650.   CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreSolidNumberOfVariables,Err)
    651. 

  651.   CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreSolidNumberOfComponents,Err)
    652. 

  652.   CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
    653. 

  653.   CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
    654. 

  654.   CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
    655. 

  655.   CALL CMISSField_CreateFinish(FibreFieldSolid,Err)
    656. 

  656. 

  657.   ! end Solid
    658. 

  658.   !--------------------------------------------------------------------------------------------------------------------------------
    659. 

  659. 

  660.   !
    661. 

  661.   !================================================================================================================================
    662. 

  662.   !
    663. 

  663. 

  664.   !EQUATIONS SETS
    665. 

  665.   DO icompartment = 1,Ncompartments
    666. 

  666.     EquationsSetFieldUserNumberDarcy = 100_CMISSIntg+icompartment
    667. 

  667.     EquationsSetUserNumberDarcy = 200_CMISSIntg+icompartment
    668. 

  668.   !Create the equations set for ALE Darcy
    669. 

  669.     CALL CMISSField_Initialise(EquationsSetFieldDarcy(icompartment),Err)
    670. 

  670.     CALL CMISSEquationsSet_Initialise(EquationsSetDarcy(icompartment),Err)
    671. 

  671.     CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricField, &
    672. 

  672.       & CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
    673. 

  673.       & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
    674. 

  674.       & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy(icompartment),EquationsSetDarcy(icompartment),Err)
    675. 

  675.     !Finish creating the equations set
    676. 

  676.     CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy(icompartment),Err)
    677. 

  677.     !Set the values for the equations set field to be the current compartment number (1 - N), and the total number of compartments (N)
    678. 

  678.     CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    679. 

  679.       & CMISS_FIELD_VALUES_SET_TYPE,1,icompartment,Err)
    680. 

  680.     CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    681. 

  681.       & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
    682. 

  682.   ENDDO
    683. 

  683. 

  684.   !--------------------------------------------------------------------------------------------------------------------------------
    685. 

  685.   ! Solid
    686. 

  686. 

  687.   !Create the equations_set
    688. 

  688.   CALL CMISSField_Initialise(EquationsSetFieldSolid,Err)
    689. 

  689.   CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err)
    690. 

  690.   CALL CMISSEquationsSet_CreateStart(EquationSetSolidUserNumber,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
    691. 

  691.     & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,&
    692. 

  692.     & EquationsSetFieldSolidUserNumber,EquationsSetFieldSolid,EquationsSetSolid,Err)
    693. 

  693.   CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err)
    694. 

  694.   !Set the values for the equations set field to be the current compartment number (O for the finite elasticity equations_set), and the total number of compartments (N)
    695. 

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

  696.   CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    697. 

  697.      & CMISS_FIELD_VALUES_SET_TYPE,1,0_CMISSIntg,Err)
    698. 

  698.   CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    699. 

  699.      & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
    700. 

  700. 

  701. 

  702.   ! end Solid
    703. 

  703.   !--------------------------------------------------------------------------------------------------------------------------------
    704. 

  704.   !--------------------------------------------------------------------------------------------------------------------------------
    705. 

  705.   ! Solid Materials Field
    706. 

  706. 

  707.   !Create a material field and attach it to the geometric field
    708. 

  708.   CALL CMISSField_Initialise(MaterialFieldSolid,Err)
    709. 

  709.   !
    710. 

  710.   CALL CMISSField_CreateStart(FieldMaterialSolidUserNumber,Region,MaterialFieldSolid,Err)
    711. 

  711.   !
    712. 

  712.   CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err)
    713. 

  713.   CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err)
    714. 

  714.   CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err)
    715. 

  715.   CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialSolidNumberOfVariables,Err)
    716. 

  716.   CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialSolidNumberOfComponents,Err)
    717. 

  717.   CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
    718. 

  718.   CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
    719. 

  719.   CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
    720. 

  720.   !
    721. 

  721.   CALL CMISSField_CreateFinish(MaterialFieldSolid,Err)
    722. 

  722. 

  723.   !Set material parameters
    724. 

  724.   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
    725. 

  725.     & 2.0_CMISSDP,Err)
    726. 

  726. !   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,2.0e3_CMISSDP,Err)
    727. 

  727.   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2, &
    728. 

  728.     & 6.0_CMISSDP,Err)
    729. 

  729. !   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,33.0_CMISSDP,Err)
    730. 

  730.   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3, &
    731. 

  731.     & 10.0_CMISSDP,Err)
    732. 

  732. 

  733.   CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialSolidUserNumber,MaterialFieldSolid,Err)
    734. 

  734.   CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err)
    735. 

  735. 

  736.   ! end Solid
    737. 

  737.   !--------------------------------------------------------------------------------------------------------------------------------
    738. 

  738. 

  739. 

  740.   !--------------------------------------------------------------------------------------------------------------------------------
    741. 

  741.   ! Solid
    742. 

  742. 

  743.   !Create a dependent field with two variables and four components
    744. 

  744.   CALL CMISSField_Initialise(DependentFieldSolid,Err)
    745. 

  745.   !
    746. 

  746.   CALL CMISSField_CreateStart(FieldDependentSolidUserNumber,Region,DependentFieldSolid,Err)
    747. 

  747.   !
    748. 

  748.   CALL CMISSField_TypeSet(DependentFieldSolid,CMISS_FIELD_GENERAL_TYPE,Err)
    749. 

  749.   CALL CMISSField_MeshDecompositionSet(DependentFieldSolid,Decomposition,Err)
    750. 

  750.   CALL CMISSField_GeometricFieldSet(DependentFieldSolid,GeometricFieldSolid,Err)
    751. 

  751.   CALL CMISSField_DependentTypeSet(DependentFieldSolid,CMISS_FIELD_DEPENDENT_TYPE,Err)
    752. 

  752.   !Create 2N+2 number of variables - 2 for solid, 2N for N Darcy compartments
    753. 

  753.   FieldDependentSolidNumberOfVariables=2*Ncompartments+2
    754. 

  754.   CALL CMISSField_NumberOfVariablesSet(DependentFieldSolid,FieldDependentSolidNumberOfVariables,Err)
    755. 

  755.   !create two variables for each compartment
    756. 

  756.   ALLOCATE(VariableTypes(2*Ncompartments+2))
    757. 

  757.   DO num_var=1,Ncompartments+1
    758. 

  758.      VariableTypes(2*num_var-1)=CMISS_FIELD_U_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
    759. 

  759.      VariableTypes(2*num_var)=CMISS_FIELD_DELUDELN_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
    760. 

  760.   ENDDO
    761. 

  761.   CALL CMISSField_VariableTypesSet(DependentFieldSolid,VariableTypes,Err) 
    762. 

  762. !   CALL CMISSField_VariableTypesSet(DependentFieldSolid,(/CMISS_FIELD_U_VARIABLE_TYPE, &
    763. 

  763. !     & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err)
    764. 

  764.     CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    765. 

  765.        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
    766. 

  766.     CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
    767. 

  767.        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
    768. 

  768.   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
    769. 

  769.   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
    770. 

  770.     & FieldDependentSolidNumberOfComponents,Err)
    771. 

  771. !   DO icompartment=3,2*Ncompartments+2
    772. 

  772. !     CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
    773. 

  773. !   ENDDO
    774. 

  774. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    775. 

  775. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    776. 

  776. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    777. 

  777. 

  778.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidDisplMeshComponentNumber,Err)
    779. 

  779.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidDisplMeshComponentNumber,Err)
    780. 

  780.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidDisplMeshComponentNumber,Err)
    781. 

  781.   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4, &
    782. 

  782.     & CMISS_FIELD_NODE_BASED_INTERPOLATION, &
    783. 

  783.     & Err)
    784. 

  784. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
    785. 

  785. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
    786. 

  786.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidLagrMultMeshComponentNumber,Err)
    787. 

  787. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
    788. 

  788. !     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    789. 

  789. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
    790. 

  790. !     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    791. 

  791. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
    792. 

  792. !     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    793. 

  793.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
    794. 

  794.     & SolidDisplMeshComponentNumber, &
    795. 

  795.     & Err)
    796. 

  796.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
    797. 

  797.     & SolidDisplMeshComponentNumber, &
    798. 

  798.     & Err)
    799. 

  799.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
    800. 

  800.     & SolidDisplMeshComponentNumber, &
    801. 

  801.     & Err)
    802. 

  802.   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
    803. 

  803.     & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    804. 

  804. !   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
    805. 

  805. !     & CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
    806. 

  806. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
    807. 

  807.   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
    808. 

  808.     & SolidLagrMultMeshComponentNumber, &
    809. 

  809.     & Err)
    810. 

  810.   !loop over the number of compartments
    811. 

  811.   DO icompartment=3,2*Ncompartments+2
    812. 

  812. !     CALL CMISSField_DimensionSet(DependentFieldSolid,VariableTypes(icompartment), &
    813. 

  813. !        & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err)
    814. 

  814.     CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err)
    815. 

  815.     DO componentnum=1,FieldDependentFluidNumberOfComponents-1
    816. 

  816.     !set dimension type
    817. 

  817. !     CALL CMISSField_DimensionSet(DependentField,VariableTypes(icompartment), &
    818. 

  818. !        & CMISS_FIELD_SCALAR_DIMENSION_TYPE,Err)
    819. 

  819.       CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, &
    820. 

  820.        & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    821. 

  821.       CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, & 
    822. 

  822.          & DarcyVelMeshComponentNumber,Err)
    823. 

  823.     ENDDO
    824. 

  824.       CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment), &
    825. 

  825.        & FieldDependentFluidNumberOfComponents, &
    826. 

  826.        & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
    827. 

  827. !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
    828. 

  828. !       & FieldDependentFluidNumberOfComponents,MESH_COMPONENT_NUMBER_PRESSURE,Err)
    829. 

  829.     CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), &
    830. 

  830.       & FieldDependentFluidNumberOfComponents,DarcyMassIncreaseMeshComponentNumber,Err)
    831. 

  831.     
    832. 

  832.   ENDDO
    833. 

  833. 

  834. !   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
    835. 

  835. !   CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
    836. 

  836. !   !For this equation type, MESH_COMPONENT_NUMBER_PRESSURE is actually the mass increase component as the pressure is taken from the solid equations
    837. 

  837. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    838. 

  838. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    839. 

  839. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    840. 

  840. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
    841. 

  841. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    842. 

  842. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    843. 

  843. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err)
    844. 

  844. !   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
    845. 

  845. 

  846.   !
    847. 

  847.   CALL CMISSField_CreateFinish(DependentFieldSolid,Err)
    848. 

  848.   !
    849. 

  849.   CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentSolidUserNumber,DependentFieldSolid,Err)
    850. 

  850.   CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err)
    851. 

  851. 

  852. !   !Initialise dependent field (solid displacement and pressure)
    853. 

  853. !   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS !+1
    854. 

  854. !     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    855. 

  855. !       & COMPONENT_NUMBER,INITIAL_FIELD_SOLID(COMPONENT_NUMBER),Err)
    856. 

  856. !   ENDDO
    857. 

  857. 

  858.   ! end Solid
    859. 

  859.   !--------------------------------------------------------------------------------------------------------------------------------
    860. 

  860. 

  861.   DO icompartment = 1,Ncompartments
    862. 

  862.     CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy(icompartment),FieldDependentSolidUserNumber,&
    863. 

  863.       & DependentFieldSolid,Err)
    864. 

  864.     CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy(icompartment),Err)
    865. 

  865.   ENDDO
    866. 

  866. 

  867. 

  868.   !Create the equations set dependent field variables for ALE Darcy
    869. 

  869. !   CALL CMISSField_Initialise(DependentFieldDarcy,Err)
    870. 

  870. !   CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,DependentFieldUserNumberDarcy, & ! ??? UserNumber ???
    871. 

  871. !   CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,FieldDependentSolidUserNumber, & ! ??? UserNumber ???
    872. 

  872. !     & DependentFieldSolid,Err)
    873. 

  873. ! !   !Set the mesh component to be used by the field components.
    874. 

  874. ! !   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
    875. 

  875. ! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
    876. 

  876. ! !       & MESH_COMPONENT_NUMBER_VELOCITY,Err)
    877. 

  877. ! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
    878. 

  878. ! !       & MESH_COMPONENT_NUMBER_VELOCITY,Err)
    879. 

  879. ! !   ENDDO
    880. 

  880. ! !   COMPONENT_NUMBER=NUMBER_OF_DIMENSIONS+1
    881. 

  881. ! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
    882. 

  882. ! !       & MESH_COMPONENT_NUMBER_PRESSURE,Err)
    883. 

  883. ! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
    884. 

  884. ! !       & MESH_COMPONENT_NUMBER_PRESSURE,Err)
    885. 

  885. !   !Finish the equations set dependent field variables
    886. 

  886. !   CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy,Err)
    887. 

  887. 

  888.   !Initialise dependent field (velocity components,pressure,mass increase)
    889. 

  889.   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS+1
    890. 

  890.     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    891. 

  891.       & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
    892. 

  892.     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    893. 

  893.       & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
    894. 

  894.   ENDDO
    895. 

  895.   !
    896. 

  896.   !================================================================================================================================
    897. 

  897.   !
    898. 

  898.   ALLOCATE(CouplingCoeffs(Ncompartments,Ncompartments))
    899. 

  899.   IF(Ncompartments==2)THEN
    900. 

  900.     CouplingCoeffs(1,1)=0.0E-01_CMISSDP
    901. 

  901. !     CouplingCoeffs(1,2)=-1.0E-04_CMISSDP
    902. 

  902. !     CouplingCoeffs(2,1)=-1.0E-04_CMISSDP
    903. 

  903.     CouplingCoeffs(1,2)=0.0E-01_CMISSDP
    904. 

  904.     CouplingCoeffs(2,1)=0.0E-01_CMISSDP
    905. 

  905.     CouplingCoeffs(2,2)=0.0E-01_CMISSDP
    906. 

  906.   ELSE IF(Ncompartments==3)THEN
    907. 

  907.     CouplingCoeffs(1,1)=1.0E-02_CMISSDP
    908. 

  908.     CouplingCoeffs(1,2)=1.0E-02_CMISSDP
    909. 

  909.     CouplingCoeffs(1,3)=0.0E-02_CMISSDP
    910. 

  910.     CouplingCoeffs(2,1)=1.0E-02_CMISSDP
    911. 

  911.     CouplingCoeffs(2,2)=2.0E-02_CMISSDP
    912. 

  912.     CouplingCoeffs(2,3)=1.0E-02_CMISSDP
    913. 

  913.     CouplingCoeffs(3,1)=0.0E-02_CMISSDP
    914. 

  914.     CouplingCoeffs(3,2)=1.0E-02_CMISSDP
    915. 

  915.     CouplingCoeffs(3,3)=1.0E-02_CMISSDP
    916. 

  916.   ELSE IF(Ncompartments==4)THEN
    917. 

  917.     CouplingCoeffs(1,1)=0.0E-02_CMISSDP
    918. 

  918.     CouplingCoeffs(1,2)=0.0E-02_CMISSDP
    919. 

  919.     CouplingCoeffs(1,3)=0.0E-02_CMISSDP
    920. 

  920.     CouplingCoeffs(1,4)=0.0E-02_CMISSDP
    921. 

  921.     CouplingCoeffs(2,1)=0.0E-02_CMISSDP
    922. 

  922.     CouplingCoeffs(2,2)=0.0E-02_CMISSDP
    923. 

  923.     CouplingCoeffs(2,3)=0.0E-02_CMISSDP
    924. 

  924.     CouplingCoeffs(2,4)=0.0E-02_CMISSDP
    925. 

  925.     CouplingCoeffs(3,1)=0.0E-02_CMISSDP
    926. 

  926.     CouplingCoeffs(3,2)=0.0E-02_CMISSDP
    927. 

  927.     CouplingCoeffs(3,3)=0.0E-02_CMISSDP
    928. 

  928.     CouplingCoeffs(3,4)=0.0E-02_CMISSDP
    929. 

  929.     CouplingCoeffs(4,1)=0.0E-02_CMISSDP
    930. 

  930.     CouplingCoeffs(4,2)=0.0E-02_CMISSDP
    931. 

  931.     CouplingCoeffs(4,3)=0.0E-02_CMISSDP
    932. 

  932.     CouplingCoeffs(4,4)=0.0E-02_CMISSDP
    933. 

  933.   ELSE IF(Ncompartments==5)THEN
    934. 

  934.     CouplingCoeffs(1,1)=0.0E-02_CMISSDP
    935. 

  935.     CouplingCoeffs(1,2)=0.0E-02_CMISSDP
    936. 

  936.     CouplingCoeffs(1,3)=0.0E-02_CMISSDP
    937. 

  937.     CouplingCoeffs(1,4)=0.0E-02_CMISSDP
    938. 

  938.     CouplingCoeffs(1,5)=0.0E-02_CMISSDP
    939. 

  939.     CouplingCoeffs(2,1)=0.0E-02_CMISSDP
    940. 

  940.     CouplingCoeffs(2,2)=0.0E-02_CMISSDP
    941. 

  941.     CouplingCoeffs(2,3)=0.0E-02_CMISSDP
    942. 

  942.     CouplingCoeffs(2,4)=0.0E-02_CMISSDP
    943. 

  943.     CouplingCoeffs(2,5)=0.0E-02_CMISSDP
    944. 

  944.     CouplingCoeffs(3,1)=0.0E-02_CMISSDP
    945. 

  945.     CouplingCoeffs(3,2)=0.0E-02_CMISSDP
    946. 

  946.     CouplingCoeffs(3,3)=0.0E-02_CMISSDP
    947. 

  947.     CouplingCoeffs(3,4)=0.0E-02_CMISSDP
    948. 

  948.     CouplingCoeffs(3,5)=0.0E-02_CMISSDP
    949. 

  949.     CouplingCoeffs(4,1)=0.0E-02_CMISSDP
    950. 

  950.     CouplingCoeffs(4,2)=0.0E-02_CMISSDP
    951. 

  951.     CouplingCoeffs(4,3)=0.0E-02_CMISSDP
    952. 

  952.     CouplingCoeffs(4,4)=0.0E-02_CMISSDP
    953. 

  953.     CouplingCoeffs(4,5)=0.0E-02_CMISSDP
    954. 

  954.     CouplingCoeffs(5,1)=0.0E-02_CMISSDP
    955. 

  955.     CouplingCoeffs(5,2)=0.0E-02_CMISSDP
    956. 

  956.     CouplingCoeffs(5,3)=0.0E-02_CMISSDP
    957. 

  957.     CouplingCoeffs(5,4)=0.0E-02_CMISSDP
    958. 

  958.     CouplingCoeffs(5,5)=0.0E-02_CMISSDP
    959. 

  959.   ELSE
    960. 

  960.     write(*,*) "Can't initialise coupling coefficients array."
    961. 

  961.   ENDIF
    962. 

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

  963.   Nparams=3
    964. 

  964.   ALLOCATE(ConstitutiveParams(Ncompartments,Nparams))
    965. 

  965.   IF(Ncompartments==2)THEN
    966. 

  966.     ConstitutiveParams(1,1)=10.0E-01_CMISSDP
    967. 

  967.     ConstitutiveParams(1,2)=10.0E-01_CMISSDP
    968. 

  968.     ConstitutiveParams(1,3)=10.0E-01_CMISSDP
    969. 

  969.     ConstitutiveParams(2,1)=10.0E-01_CMISSDP
    970. 

  970.     ConstitutiveParams(2,2)=10.0E-01_CMISSDP
    971. 

  971.     ConstitutiveParams(2,3)=10.0E-01_CMISSDP
    972. 

  972.   ELSE IF(Ncompartments==3)THEN
    973. 

  973.     ConstitutiveParams(1,1)=1.0E-02_CMISSDP
    974. 

  974.     ConstitutiveParams(1,2)=1.0E-02_CMISSDP
    975. 

  975.     ConstitutiveParams(1,3)=0.0E-02_CMISSDP
    976. 

  976.     ConstitutiveParams(2,1)=1.0E-02_CMISSDP
    977. 

  977.     ConstitutiveParams(2,2)=2.0E-02_CMISSDP
    978. 

  978.     ConstitutiveParams(2,3)=1.0E-02_CMISSDP
    979. 

  979.     ConstitutiveParams(3,1)=0.0E-02_CMISSDP
    980. 

  980.     ConstitutiveParams(3,2)=1.0E-02_CMISSDP
    981. 

  981.     ConstitutiveParams(3,3)=1.0E-02_CMISSDP
    982. 

  982.   ELSE IF(Ncompartments==4)THEN
    983. 

  983.     ConstitutiveParams(1,1)=0.0E-02_CMISSDP
    984. 

  984.     ConstitutiveParams(1,2)=0.0E-02_CMISSDP
    985. 

  985.     ConstitutiveParams(1,3)=0.0E-02_CMISSDP
    986. 

  986.     ConstitutiveP…
    987.