PageRenderTime 60ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 0ms

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

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

  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.     ConstitutiveParams(2,1)=0.0E-02_CMISSDP
    987. 

  987.     ConstitutiveParams(2,2)=0.0E-02_CMISSDP
    988. 

  988.     ConstitutiveParams(2,3)=0.0E-02_CMISSDP
    989. 

  989.     ConstitutiveParams(3,1)=0.0E-02_CMISSDP
    990. 

  990.     ConstitutiveParams(3,2)=0.0E-02_CMISSDP
    991. 

  991.     ConstitutiveParams(3,3)=0.0E-02_CMISSDP
    992. 

  992.     ConstitutiveParams(4,1)=0.0E-02_CMISSDP
    993. 

  993.     ConstitutiveParams(4,2)=0.0E-02_CMISSDP
    994. 

  994.     ConstitutiveParams(4,3)=0.0E-02_CMISSDP
    995. 

  995.   ELSE IF(Ncompartments==5)THEN
    996. 

  996.     ConstitutiveParams(1,1)=0.0E-02_CMISSDP
    997. 

  997.     ConstitutiveParams(1,2)=0.0E-02_CMISSDP
    998. 

  998.     ConstitutiveParams(1,3)=0.0E-02_CMISSDP
    999. 

  999.     ConstitutiveParams(2,1)=0.0E-02_CMISSDP
    1000. 

  1000.     ConstitutiveParams(2,2)=0.0E-02_CMISSDP
    1001. 

  1001.     ConstitutiveParams(2,3)=0.0E-02_CMISSDP
    1002. 

  1002.     ConstitutiveParams(3,1)=0.0E-02_CMISSDP
    1003. 

  1003.     ConstitutiveParams(3,2)=0.0E-02_CMISSDP
    1004. 

  1004.     ConstitutiveParams(3,3)=0.0E-02_CMISSDP
    1005. 

  1005.     ConstitutiveParams(4,1)=0.0E-02_CMISSDP
    1006. 

  1006.     ConstitutiveParams(4,2)=0.0E-02_CMISSDP
    1007. 

  1007.     ConstitutiveParams(4,3)=0.0E-02_CMISSDP
    1008. 

  1008.     ConstitutiveParams(5,1)=0.0E-02_CMISSDP
    1009. 

  1009.     ConstitutiveParams(5,2)=0.0E-02_CMISSDP
    1010. 

  1010.     ConstitutiveParams(5,3)=0.0E-02_CMISSDP
    1011. 

  1011.   ELSE
    1012. 

  1012.     write(*,*) "Can't initialise constitutive parameters array."
    1013. 

  1013.   ENDIF
    1014. 

  1014.   !
    1015. 

  1015.   !================================================================================================================================
    1016. 

  1016.   !
    1017. 

  1017. 

  1018.   !MATERIALS FIELDS
    1019. 

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

  1020.   !It also contains a V variable type to store the coupling coefficients 
    1021. 

  1021.   DO icompartment = 1,Ncompartments
    1022. 

  1022.     MaterialsFieldUserNumberDarcy = 400+icompartment
    1023. 

  1023.     CALL CMISSField_Initialise(MaterialsFieldDarcy(icompartment),Err)
    1024. 

  1024.     CALL CMISSField_CreateStart(MaterialsFieldUserNumberDarcy,Region,MaterialsFieldDarcy(icompartment),Err)
    1025. 

  1025.     CALL CMISSField_TypeSet(MaterialsFieldDarcy(icompartment),CMISS_FIELD_MATERIAL_TYPE,Err)
    1026. 

  1026.     CALL CMISSField_MeshDecompositionSet(MaterialsFieldDarcy(icompartment),Decomposition,Err)
    1027. 

  1027.     CALL CMISSField_GeometricFieldSet(MaterialsFieldDarcy(icompartment),GeometricField,Err)
    1028. 

  1028.     CALL CMISSField_DependentTypeSet(MaterialsFieldDarcy(icompartment),CMISS_FIELD_INDEPENDENT_TYPE,Err)
    1029. 

  1029.     CALL CMISSField_NumberOfVariablesSet(MaterialsFieldDarcy(icompartment),3,Err)
    1030. 

  1030.     CALL CMISSField_VariableTypesSet(MaterialsFieldDarcy(icompartment),[CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,&
    1031. 

  1031.        & CMISS_FIELD_U1_VARIABLE_TYPE],Err) 
    1032. 

  1032.     CALL CMISSField_NumberOfComponentsSet(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE,7,Err)
    1033. 

  1033.     CALL CMISSField_NumberOfComponentsSet(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE,Ncompartments,Err)
    1034. 

  1034.     CALL CMISSField_NumberOfComponentsSet(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U1_VARIABLE_TYPE,Nparams,Err)
    1035. 

  1035. 

  1036.     CALL CMISSField_CreateFinish(MaterialsFieldDarcy(icompartment),Err)
    1037. 

  1037.     CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy(icompartment),MaterialsFieldUserNumberDarcy,&
    1038. 

  1038.          & MaterialsFieldDarcy(icompartment),Err)
    1039. 

  1039.     CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy(icompartment),Err)
    1040. 

  1040.   END DO
    1041. 

  1041.   DO icompartment = 1,Ncompartments
    1042. 

  1042.     CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    1043. 

  1043.       & CMISS_FIELD_VALUES_SET_TYPE, &
    1044. 

  1044.       & MaterialsFieldUserNumberDarcyPorosity,POROSITY_PARAM_DARCY,Err)
    1045. 

  1045.     CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    1046. 

  1046.       & CMISS_FIELD_VALUES_SET_TYPE, &
    1047. 

  1047.       & MaterialsFieldUserNumberDarcyPermOverVis,PERM_OVER_VIS_PARAM_DARCY,Err)
    1048. 

  1048.   END DO
    1049. 

  1049.   DO icompartment = 1, Ncompartments
    1050. 

  1050.     DO COMPONENT_NUMBER=1, Ncompartments
    1051. 

  1051.       CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
    1052. 

  1052.          & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
    1053. 

  1053. !         CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
    1054. 

  1054. !           & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
    1055. 

  1055.     END DO
    1056. 

  1056.   END DO
    1057. 

  1057.   DO icompartment = 1, Ncompartments
    1058. 

  1058.     DO COMPONENT_NUMBER=1,Nparams
    1059. 

  1059.       CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U1_VARIABLE_TYPE, &
    1060. 

  1060.          & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,ConstitutiveParams(icompartment,COMPONENT_NUMBER),Err)
    1061. 

  1061. !         CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
    1062. 

  1062. !           & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
    1063. 

  1063.     END DO
    1064. 

  1064.   END DO
    1065. 

  1065.   !Create the equations set materials field variables for ALE Darcy
    1066. 

  1066. !   CALL CMISSField_Initialise(MaterialsFieldDarcy,Err)
    1067. 

  1067. !   CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy,MaterialsFieldUserNumberDarcy, &
    1068. 

  1068. !     & MaterialsFieldDarcy,Err)
    1069. 

  1069. !   !Finish the equations set materials field variables
    1070. 

  1070. !   CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy,Err)
    1071. 

  1071. !   CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    1072. 

  1072. !     & MaterialsFieldUserNumberDarcyPorosity,POROSITY_PARAM_DARCY,Err)
    1073. 

  1073. !   CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
    1074. 

  1074. !     & MaterialsFieldUserNumberDarcyPermOverVis,PERM_OVER_VIS_PARAM_DARCY,Err)
    1075. 

  1075.   !Create the equations set materials field variables for deformation-dependent material properties
    1076. 

  1076. 

  1077.   !
    1078. 

  1078.   !================================================================================================================================
    1079. 

  1079.   !
    1080. 

  1080.   DO icompartment=1,Ncompartments
    1081. 

  1081.      SourceFieldDarcyUserNumber = 450+icompartment
    1082. 

  1082.     CALL CMISSField_Initialise(SourceFieldDarcy(icompartment),Err)
    1083. 

  1083.     CALL CMISSEquationsSet_SourceCreateStart(EquationsSetDarcy(icompartment),SourceFieldDarcyUserNumber,&
    1084. 

  1084.         & SourceFieldDarcy(icompartment),Err)
    1085. 

  1085.     CALL CMISSEquationsSet_SourceCreateFinish(EquationsSetDarcy(icompartment),Err)
    1086. 

  1086.   ENDDO
    1087. 

  1087.   !EQUATIONS
    1088. 

  1088. 

  1089.   !INDEPENDENT FIELD Darcy for storing BC flags
    1090. 

  1090.   DO icompartment=1,Ncompartments
    1091. 

  1091.      IndependentFieldDarcyUserNumber = 470+icompartment
    1092. 

  1092.     CALL CMISSField_Initialise(IndependentFieldDarcy(icompartment),Err)
    1093. 

  1093.     CALL CMISSEquationsSet_IndependentCreateStart(EquationsSetDarcy(icompartment),IndependentFieldDarcyUserNumber, &
    1094. 

  1094.       & IndependentFieldDarcy(icompartment),Err)
    1095. 

  1095. 

  1096.     CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE,1, &
    1097. 

  1097.       & DarcyVelMeshComponentNumber,Err)
    1098. 

  1098.     CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE,2, &
    1099. 

  1099.       & DarcyVelMeshComponentNumber,Err)
    1100. 

  1100.     CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE,3, &
    1101. 

  1101.       & DarcyVelMeshComponentNumber,Err)
    1102. 

  1102. ! !   CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,4,DarcyMassIncreaseMeshComponentNumber,Err)
    1103. 

  1103. 

  1104.     CALL CMISSEquationsSet_IndependentCreateFinish(EquationsSetDarcy(icompartment),Err)
    1105. 

  1105. 

  1106.     CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    1107. 

  1107.       & CMISS_FIELD_VALUES_SET_TYPE,1,0.0_CMISSDP,Err)
    1108. 

  1108.     CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
    1109. 

  1109.       & CMISS_FIELD_VALUES_SET_TYPE,2,0.0_CMISSDP,Err)
    1110. 

  1110.     CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, & 
    1111. 

  1111.       & CMISS_FIELD_VALUES_SET_TYPE,3,0.0_CMISSDP,Err)
    1112. 

  1112.   ENDDO
    1113. 

  1113.  
    1114. 

  1114.   DO icompartment=1,Ncompartments
    1115. 

  1115.     !Create the equations set equations
    1116. 

  1116.     CALL CMISSEquations_Initialise(EquationsDarcy(icompartment),Err)
    1117. 

  1117.     CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy(icompartment),EquationsDarcy(icompartment),Err)
    1118. 

  1118.     !Set the equations matrices sparsity type
    1119. 

  1119.     CALL CMISSEquations_SparsityTypeSet(EquationsDarcy(icompartment),CMISS_EQUATIONS_SPARSE_MATRICES,Err)
    1120. 

  1120.   !   !Set the equations lumping type
    1121. 

  1121.   !   CALL CMISSEquations_LumpingTypeSet(EquationsDarcy,CMISS_EQUATIONS_UNLUMPED_MATRICES,Err)
    1122. 

  1122.     !Set the equations set output
    1123. 

  1123.     CALL CMISSEquations_OutputTypeSet(EquationsDarcy(icompartment),EQUATIONS_DARCY_OUTPUT,Err)
    1124. 

  1124.   !Finish the equations set equations
    1125. 

  1125.     CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetDarcy(icompartment),Err)
    1126. 

  1126.   ENDDO
    1127. 

  1127.   !--------------------------------------------------------------------------------------------------------------------------------
    1128. 

  1128.   ! Solid
    1129. 

  1129. 

  1130.   !Create the equations set equations
    1131. 

  1131.   CALL CMISSEquations_Initialise(EquationsSolid,Err)
    1132. 

  1132.   CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetSolid,EquationsSolid,Err)
    1133. 

  1133.   CALL CMISSEquations_SparsityTypeSet(EquationsSolid,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
    1134. 

  1134.   CALL CMISSEquations_OutputTypeSet(EquationsSolid,CMISS_EQUATIONS_NO_OUTPUT,Err)
    1135. 

  1135.   CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetSolid,Err)
    1136. 

  1136. 

  1137.   ! end Solid
    1138. 

  1138.   !--------------------------------------------------------------------------------------------------------------------------------
    1139. 

  1139. 

  1140.   !
    1141. 

  1141.   !================================================================================================================================
    1142. 

  1142.   !
    1143. 

  1143. 

  1144.   !--------------------------------------------------------------------------------------------------------------------------------
    1145. 

  1145.   ! Solid
    1146. 

  1146. 

  1147.   !Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure
    1148. 

  1148.   CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    1149. 

  1149.     & CMISS_FIELD_VALUES_SET_TYPE, &
    1150. 

  1150.     & 1,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,Err)
    1151. 

  1151.   CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    1152. 

  1152.     & CMISS_FIELD_VALUES_SET_TYPE, &
    1153. 

  1153.     & 2,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,Err)
    1154. 

  1154.   CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
    1155. 

  1155.     & CMISS_FIELD_VALUES_SET_TYPE, &
    1156. 

  1156.     & 3,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,Err)
    1157. 

  1157.   CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4, &
    1158. 

  1158.     & 0.0_CMISSDP, &
    1159. 

  1159.     & Err)
    1160. 

  1160. 

  1161.   ! end Solid
    1162. 

  1162.   !--------------------------------------------------------------------------------------------------------------------------------
    1163. 

  1163. 

  1164.   !
    1165. 

  1165.   !================================================================================================================================
    1166. 

  1166.   !
    1167. 

  1167. 

  1168.   !PROBLEMS
    1169. 

  1169. 

  1170.   CALL CMISSProblem_Initialise(Problem,Err)
    1171. 

  1171.   CALL CMISSControlLoop_Initialise(ControlLoop,Err)
    1172. 

  1172.   CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
    1173. 

  1173.   CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_MULTI_PHYSICS_CLASS,CMISS_PROBLEM_FINITE_ELASTICITY_DARCY_TYPE, &
    1174. 

  1174.     & CMISS_PROBLEM_QUASISTATIC_ELASTICITY_TRANSIENT_DARCY_SUBTYPE,Err)
    1175. 

  1175.   CALL CMISSProblem_CreateFinish(Problem,Err)
    1176. 

  1176. 

  1177.   CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
    1178. 

  1178.   CALL CMISSProblem_ControlLoopGet(Problem,CMISS_CONTROL_LOOP_NODE,ControlLoop,Err)
    1179. 

  1179. !   CALL CMISSControlLoop_MaximumIterationsSet(ControlLoop,2,Err)
    1180. 

  1180.   CALL CMISSControlLoop_TimesSet(ControlLoop,DYNAMIC_SOLVER_DARCY_START_TIME,DYNAMIC_SOLVER_DARCY_STOP_TIME, &
    1181. 

  1181.     & DYNAMIC_SOLVER_DARCY_TIME_INCREMENT,Err)
    1182. 

  1182.   CALL CMISSControlLoop_TimeOutputSet(ControlLoop,DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY,Err)
    1183. 

  1183. !   CALL CMISSControlLoop_OutputTypeSet(ControlLoop,CMISS_CONTROL_LOOP_PROGRESS_OUTPUT,Err)
    1184. 

  1184.   CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
    1185. 

  1185. 

  1186. 

  1187. 

  1188.   !
    1189. 

  1189.   !================================================================================================================================
    1190. 

  1190.   !
    1191. 

  1191. 

  1192.   !SOLVERS
    1193. 

  1193. 

  1194.   !Start the creation of the problem solvers
    1195. 

  1195.   CALL CMISSSolver_Initialise(SolverSolid,Err)
    1196. 

  1196.   CALL CMISSSolver_Initialise(DynamicSolverDarcy,Err)
    1197. 

  1197.   CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
    1198. 

  1198. 

  1199.   CALL CMISSProblem_SolversCreateStart(Problem,Err)
    1200. 

  1200. 

  1201.   !--------------------------------------------------------------------------------------------------------------------------------
    1202. 

  1202.   ! Solid
    1203. 

  1203. 

  1204.   !Get the finite elasticity solver
    1205. 

  1205.   CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
    1206. 

  1206.     & SolverSolidIndex,SolverSolid,Err)
    1207. 

  1207.   CALL CMISSSolver_OutputTypeSet(SolverSolid,CMISS_SOLVER_PROGRESS_OUTPUT,Err)
    1208. 

  1208. !   CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_FD_CALCULATED,Err)
    1209. 

  1209.   CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_EQUATIONS_CALCULATED,Err)
    1210. 

  1210. 

  1211.   CALL CMISSSolver_NewtonAbsoluteToleranceSet(SolverSolid,ABSOLUTE_TOLERANCE,Err)
    1212. 

  1212.   CALL CMISSSolver_NewtonRelativeToleranceSet(SolverSolid,RELATIVE_TOLERANCE,Err)
    1213. 

  1213.   CALL CMISSSolver_NewtonMaximumIterationsSet(SolverSolid,MAXIMUM_ITERATIONS,Err)
    1214. 

  1214. 

  1215. !   CALL CMISSSolverNonLinearTypeSet(SolverSolid,CMISS_SOLVER_NONLINEAR_NEWTON,Err)
    1216. 

  1216. !   CALL CMISSSolver_LibraryTypeSet(SolverSolid,CMISS_SOLVER_PETSC_LIBRARY,Err)
    1217. 

  1217. 

  1218. !   CALL CMISSSolver_NewtonLinearSolverGet(SolverSolid,LinearSolverSolid,Err)
    1219. 

  1219. !   CALL CMISSSolver_LinearTypeSet(LinearSolverSolid,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
    1220. 

  1220. 

  1221.   ! end Solid
    1222. 

  1222.   !--------------------------------------------------------------------------------------------------------------------------------
    1223. 

  1223. 

  1224.   !Get the Darcy solver
    1225. 

  1225.   CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
    1226. 

  1226.     & SolverDarcyIndex,DynamicSolverDarcy,Err)
    1227. 

  1227.   !Set the output type
    1228. 

  1228.   CALL CMISSSolver_OutputTypeSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE,Err)
    1229. 

  1229.   !Set theta
    1230. 

  1230.   CALL CMISSSolver_DynamicThetaSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_THETA,Err)
    1231. 

  1231.   !CALL CMISSSolverDynamicDynamicSet(DynamicSolverDarcy,.TRUE.,Err)
    1232. 

  1232.   !Get the dynamic linear solver
    1233. 

  1233.   CALL CMISSSolver_DynamicLinearSolverGet(DynamicSolverDarcy,LinearSolverDarcy,Err)
    1234. 

  1234.   !Set the solver settings
    1235. 

  1235.   IF(LINEAR_SOLVER_DARCY_DIRECT_FLAG) THEN
    1236. 

  1236.     CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
    1237. 

  1237.     CALL CMISSSolver_LibraryTypeSet(LinearSolverDarcy,CMISS_SOLVER_MUMPS_LIBRARY,Err)
    1238. 

  1238.   ELSE
    1239. 

  1239.     CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
    1240. 

  1240.     CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverDarcy,MAXIMUM_ITERATIONS,Err)
    1241. 

  1241.     CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverDarcy,DIVERGENCE_TOLERANCE,Err)
    1242. 

  1242.     CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverDarcy,RELATIVE_TOLERANCE,Err)
    1243. 

  1243.     CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverDarcy,ABSOLUTE_TOLERANCE,Err)
    1244. 

  1244.     CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverDarcy,RESTART_VALUE,Err)
    1245. 

  1245.   ENDIF
    1246. 

  1246. 

  1247.   !Finish the creation of the problem solver
    1248. 

  1248.   CALL CMISSProblem_SolversCreateFinish(Problem,Err)
    1249. 

  1249. 

  1250.   !
    1251. 

  1251.   !================================================================================================================================
    1252. 

  1252.   !
    1253. 

  1253. 

  1254.   !SOLVER EQUATIONS
    1255. 

  1255. 

  1256.   !Start the creation of the problem solver equations
    1257. 

  1257.   CALL CMISSSolver_Initialise(SolverSolid,Err)
    1258. 

  1258.   CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
    1259. 

  1259. 

  1260.   CALL CMISSSolverEquations_Initialise(SolverEquationsSolid,Err)
    1261. 

  1261.   CALL CMISSSolverEquations_Initialise(SolverEquationsDarcy,Err)
    1262. 

  1262. 

  1263.   CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
    1264. 

  1264.   !
    1265. 

  1265.   !Get the finite elasticity solver equations
    1266. 

  1266.   CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
    1267. 

  1267.     & SolverSolidIndex,SolverSolid,Err)
    1268. 

  1268.   CALL CMISSSolver_SolverEquationsGet(SolverSolid,SolverEquationsSolid,Err)
    1269. 

  1269.   CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsSolid,CMISS_SOLVER_SPARSE_MATRICES,Err)
    1270. 

  1270.   CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsSolid,EquationsSetSolid,EquationsSetIndex,Err)
    1271. 

  1271.   !
    1272. 

  1272.   !Get the Darcy solver equations
    1273. 

  1273.   CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
    1274. 

  1274.     & SolverDarcyIndex,LinearSolverDarcy,Err)
    1275. 

  1275.   CALL CMISSSolver_SolverEquationsGet(LinearSolverDarcy,SolverEquationsDarcy,Err)
    1276. 

  1276.   CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsDarcy,CMISS_SOLVER_SPARSE_MATRICES,Err)
    1277. 

  1277.   DO icompartment=1,Ncompartments
    1278. 

  1278.     CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsDarcy,EquationsSetDarcy(icompartment),EquationsSetIndex,Err)
    1279. 

  1279.   ENDDO
    1280. 

  1280.   !
    1281. 

  1281.   !Finish the creation of the problem solver equations
    1282. 

  1282.   CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
    1283. 

  1283. !   ENDDO
    1284. 

  1284. 

  1285. 

  1286. 

  1287. 

  1288.   !Prescribe boundary conditions (absolute nodal parameters)
    1289. 

  1289.   !Solid is computed in absolute position, rather than displacement. Thus BCs for absolute position
    1290. 

  1290.    ! ASSIGN BOUNDARY CONDITIONS - SOLID (absolute nodal parameters)
    1291. 

  1291.   !Solid is computed in absolute position, rather than displacement. Thus BCs for absolute position
    1292. 

  1292.   CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsSolid,Err)
    1293. 

  1293.   CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsSolid,BoundaryConditionsSolid,Err)
    1294. 

  1294. 

  1295.   !Get surfaces
    1296. 

  1296.   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_FRONT_SURFACE, &
    1297. 

  1297.     & Face1Nodes,FaceXi(1),Err)
    1298. 

  1298.   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BACK_SURFACE, &
    1299. 

  1299.     & Face2Nodes,FaceXi(2),Err)
    1300. 

  1300.   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_RIGHT_SURFACE, &
    1301. 

  1301.     & Face3Nodes,FaceXi(3),Err)
    1302. 

  1302.   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_LEFT_SURFACE, &
    1303. 

  1303.     & Face4Nodes,FaceXi(4),Err)
    1304. 

  1304.   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_TOP_SURFACE, &
    1305. 

  1305.     & Face5Nodes,FaceXi(5),Err)
    1306. 

  1306.   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BOTTOM_SURFACE, &
    1307. 

  1307.     & Face6Nodes,FaceXi(6),Err)
    1308. 

  1308. 

  1309.   ! Fix the bottom in z direction
    1310. 

  1310.   DO NN=1,SIZE(Face6Nodes,1)
    1311. 

  1311.     NODE=Face6Nodes(NN)
    1312. 

  1312.     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1313. 

  1313.     IF(NodeDomain==ComputationalNodeNumber) THEN
    1314. 

  1314.       CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, &
    1315. 

  1315.         & ZCoord,Err)
    1316. 

  1316.       CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, &
    1317. 

  1317.         & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err)
    1318. 

  1318.       WRITE(*,*) "FIXING NODE",NODE,"AT BOTTOM IN Z DIRECTION"
    1319. 

  1319.     ENDIF
    1320. 

  1320.   ENDDO
    1321. 

  1321. 

  1322.   ! Fix the top in z direction
    1323. 

  1323.   DO NN=1,SIZE(Face5Nodes,1)
    1324. 

  1324.     NODE=Face5Nodes(NN)
    1325. 

  1325.     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1326. 

  1326.     IF(NodeDomain==ComputationalNodeNumber) THEN
    1327. 

  1327.       CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, &
    1328. 

  1328.         & ZCoord,Err)
    1329. 

  1329.       CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, &
    1330. 

  1330.         & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err)
    1331. 

  1331.       WRITE(*,*) "FIXING NODE",NODE,"AT TOP IN Z DIRECTION"
    1332. 

  1332.     ENDIF
    1333. 

  1333.   ENDDO
    1334. 

  1334. 

  1335.   !Fix more nodes at the bottom to stop free body motion
    1336. 

  1336.   X_FIXED=.FALSE.
    1337. 

  1337.   Y_FIXED=.FALSE.
    1338. 

  1338.   DO NN=1,SIZE(Face6Nodes,1)
    1339. 

  1339.     NODE=Face6Nodes(NN)
    1340. 

  1340.     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1341. 

  1341.     IF(NodeDomain==ComputationalNodeNumber) THEN
    1342. 

  1342.       CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,1, &
    1343. 

  1343.         & XCoord,Err)
    1344. 

  1344.       CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,2, &
    1345. 

  1345.         & YCoord,Err)
    1346. 

  1346. 

  1347.       !Fix Origin displacement in x and y (z already fixed)
    1348. 

  1348.       IF(ABS(XCoord)<1.0E-6_CMISSDP.AND.ABS(YCoord)<1.0E-6_CMISSDP) THEN
    1349. 

  1349.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,1, &
    1350. 

  1350.           & CMISS_BOUNDARY_CONDITION_FIXED,XCoord,Err)
    1351. 

  1351.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,2, &
    1352. 

  1352.           & CMISS_BOUNDARY_CONDITION_FIXED,YCoord,Err)
    1353. 

  1353.         WRITE(*,*) "FIXING ORIGIN NODE",NODE,"IN X AND Y DIRECTION"
    1354. 

  1354.         X_FIXED=.TRUE.
    1355. 

  1355.         Y_FIXED=.TRUE.
    1356. 

  1356.       ENDIF
    1357. 

  1357. 

  1358.       !Fix nodal displacements at (X_DIM,0) in y
    1359. 

  1359.       IF(ABS(XCoord - X_DIM)<1.0E-6_CMISSDP .AND. ABS(YCoord)<1.0E-6_CMISSDP) THEN
    1360. 

  1360.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,2, &
    1361. 

  1361.           & CMISS_BOUNDARY_CONDITION_FIXED,YCoord,Err)
    1362. 

  1362.         WRITE(*,*) "FIXING NODES",NODE,"AT (X_DIM,0) IN Y DIRECTION"
    1363. 

  1363.         Y_FIXED=.TRUE.
    1364. 

  1364.       ENDIF
    1365. 

  1365. 

  1366.       !Fix nodal displacements at (0,Y_DIM) in x
    1367. 

  1367.       IF(ABS(XCoord)<1.0E-6_CMISSDP .AND. ABS(YCoord - Y_DIM)<1.0E-6_CMISSDP) THEN
    1368. 

  1368.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,1, &
    1369. 

  1369.           & CMISS_BOUNDARY_CONDITION_FIXED,XCoord,Err)
    1370. 

  1370.         WRITE(*,*) "FIXING NODES",NODE,"AT (0,Y_DIM) IN X DIRECTION"
    1371. 

  1371.         X_FIXED=.TRUE.
    1372. 

  1372.       ENDIF
    1373. 

  1373. 

  1374.     ENDIF
    1375. 

  1375.   ENDDO
    1376. 

  1376. !   CALL MPI_REDUCE(X_FIXED,X_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
    1377. 

  1377. !   CALL MPI_REDUCE(Y_FIXED,Y_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
    1378. 

  1378. !   IF(ComputationalNodeNumber==0) THEN
    1379. 

  1379. !     IF(.NOT.(X_OKAY.AND.Y_OKAY)) THEN
    1380. 

  1380. !       WRITE(*,*) "Free body motion could not be prevented!"
    1381. 

  1381. !       CALL CMISSFinalise(Err)
    1382. 

  1382. !       STOP
    1383. 

  1383. !     ENDIF
    1384. 

  1384. !   ENDIF
    1385. 

  1385. 

  1386.   CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsSolid,Err)
    1387. 

  1387. 

  1388. 

  1389.   !
    1390. 

  1390.   !------------------------------------
    1391. 

  1391.   ! ASSIGN BOUNDARY CONDITIONS - FLUID
    1392. 

  1392.     !Get surfaces
    1393. 

  1393.     CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_FRONT_SURFACE, &
    1394. 

  1394.       & Face7Nodes,FaceXi(1),Err)
    1395. 

  1395.     CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BACK_SURFACE, &
    1396. 

  1396.       & Face8Nodes,FaceXi(2),Err)
    1397. 

  1397.     CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_RIGHT_SURFACE, &
    1398. 

  1398.       & Face9Nodes,FaceXi(3),Err)
    1399. 

  1399.     CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_LEFT_SURFACE, &
    1400. 

  1400.       & Face10Nodes,FaceXi(4),Err)
    1401. 

  1401.     CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_TOP_SURFACE, &
    1402. 

  1402.       & Face11Nodes,FaceXi(5),Err)
    1403. 

  1403.     CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BOTTOM_SURFACE, &
    1404. 

  1404.       & Face12Nodes,FaceXi(6),Err)
    1405. 

  1405. 

  1406.   CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsDarcy,Err)
    1407. 

  1407.   CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsDarcy,BoundaryConditionsDarcy,Err)
    1408. 

  1408.   DO icompartment=1,Ncompartments
    1409. 

  1409. 

  1410.     IF(icompartment==1) THEN
    1411. 

  1411.     ! At the top impose Darcy velocity in z direction
    1412. 

  1412.     DO NN=1,SIZE(Face11Nodes,1)
    1413. 

  1413.       NODE=Face11Nodes(NN)
    1414. 

  1414. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1415. 

  1415. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1416. 

  1416.         VALUE = -1.0_CMISSDP
    1417. 

  1417.         COMPONENT_NUMBER = 3
    1418. 

  1418.         write(*,*)'Marker 0'
    1419. 

  1419.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,1,NODE, &
    1420. 

  1420.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1421. 

  1421.         WRITE(*,*) "SPECIFIED INFLOW AT NODE",NODE,"IN Z DIRECTION"
    1422. 

  1422. 

  1423. ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,1,XCoord,Err)
    1424. 

  1424. ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,2,YCoord,Err)
    1425. 

  1425. ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,3,ZCoord,Err)
    1426. 

  1426. ! !       WRITE(*,*) "XCoord, YCoord, ZCoord = ",XCoord, YCoord, ZCoord
    1427. 

  1427. ! !     ENDIF
    1428. 

  1428.     ENDDO
    1429. 

  1429. 

  1430.     !All other faces are impermeable
    1431. 

  1431.     DO NN=1,SIZE(Face7Nodes,1)
    1432. 

  1432.       NODE=Face7Nodes(NN)
    1433. 

  1433. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1434. 

  1434. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1435. 

  1435.         VALUE = 0.0_CMISSDP
    1436. 

  1436.         COMPONENT_NUMBER = 1
    1437. 

  1437.         write(*,*)'Marker 1'
    1438. 

  1438.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,1,NODE, &
    1439. 

  1439.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1440. 

  1440.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN X DIRECTION"
    1441. 

  1441. ! !     ENDIF
    1442. 

  1442.     ENDDO
    1443. 

  1443. 

  1444.     DO NN=1,SIZE(Face8Nodes,1)
    1445. 

  1445.       NODE=Face8Nodes(NN)
    1446. 

  1446. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1447. 

  1447. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1448. 

  1448.         VALUE = 0.0_CMISSDP
    1449. 

  1449.         COMPONENT_NUMBER = 1
    1450. 

  1450.         write(*,*)'Marker 2'
    1451. 

  1451.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,1,NODE, &
    1452. 

  1452.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1453. 

  1453.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN X DIRECTION"
    1454. 

  1454. ! !     ENDIF
    1455. 

  1455.     ENDDO
    1456. 

  1456. 

  1457.     DO NN=1,SIZE(Face9Nodes,1)
    1458. 

  1458.       NODE=Face9Nodes(NN)
    1459. 

  1459. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1460. 

  1460. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1461. 

  1461.         VALUE = 0.0_CMISSDP
    1462. 

  1462.         COMPONENT_NUMBER = 2
    1463. 

  1463.         write(*,*)'Marker 3'
    1464. 

  1464.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,1,NODE, &
    1465. 

  1465.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1466. 

  1466.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN Y DIRECTION"
    1467. 

  1467. ! !     ENDIF
    1468. 

  1468.     ENDDO
    1469. 

  1469. 

  1470.     DO NN=1,SIZE(Face10Nodes,1)
    1471. 

  1471.       NODE=Face10Nodes(NN)
    1472. 

  1472. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1473. 

  1473. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1474. 

  1474.         VALUE = 0.0_CMISSDP
    1475. 

  1475.         COMPONENT_NUMBER = 2
    1476. 

  1476.         write(*,*)'Marker 4'
    1477. 

  1477.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,1,NODE, &
    1478. 

  1478.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1479. 

  1479.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN Y DIRECTION"
    1480. 

  1480. ! !     ENDIF
    1481. 

  1481.     ENDDO
    1482. 

  1482. 

  1483.     DO NN=1,SIZE(Face12Nodes,1)
    1484. 

  1484.       NODE=Face12Nodes(NN)
    1485. 

  1485. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1486. 

  1486. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1487. 

  1487.         VALUE = 0.0_CMISSDP
    1488. 

  1488.         COMPONENT_NUMBER = 3
    1489. 

  1489.         write(*,*)'Marker 5'
    1490. 

  1490.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,1,NODE, &
    1491. 

  1491.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1492. 

  1492.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN Z DIRECTION"
    1493. 

  1493. ! !     ENDIF
    1494. 

  1494.     ENDDO
    1495. 

  1495.     ENDIF
    1496. 

  1496.     !Compartment TWO boundary conditions!
    1497. 

  1497.     ! At the top impose Darcy velocity in z direction
    1498. 

  1498.     IF(icompartment==2) THEN
    1499. 

  1499.     DO NN=1,SIZE(Face11Nodes,1)
    1500. 

  1500.       NODE=Face11Nodes(NN)
    1501. 

  1501. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1502. 

  1502. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1503. 

  1503.         VALUE = -1.0_CMISSDP
    1504. 

  1504.         COMPONENT_NUMBER = 3
    1505. 

  1505.         write(*,*)'Marker 0'
    1506. 

  1506.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,1,1,NODE, &
    1507. 

  1507.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1508. 

  1508.         WRITE(*,*) "SPECIFIED INFLOW AT NODE",NODE,"IN Z DIRECTION"
    1509. 

  1509. 

  1510. ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,1,XCoord,Err)
    1511. 

  1511. ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,2,YCoord,Err)
    1512. 

  1512. ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,3,ZCoord,Err)
    1513. 

  1513. ! !       WRITE(*,*) "XCoord, YCoord, ZCoord = ",XCoord, YCoord, ZCoord
    1514. 

  1514. ! !     ENDIF
    1515. 

  1515.     ENDDO
    1516. 

  1516. !     DO NN=1,SIZE(Face11Nodes,1)
    1517. 

  1517. !       NODE=Face11Nodes(NN)
    1518. 

  1518. ! ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1519. 

  1519. ! ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1520. 

  1520. !         VALUE = 0.0_CMISSDP
    1521. 

  1521. !         COMPONENT_NUMBER = 3
    1522. 

  1522. !         write(*,*)'Marker 0'
    1523. 

  1523. !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,1,1,NODE,COMPONENT_NUMBER,&
    1524. 

  1524. !           & CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1525. 

  1525. !         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN Z DIRECTION"
    1526. 

  1526. !
    1527. 

  1527. ! ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,1,XCoord,Err)
    1528. 

  1528. ! ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,2,YCoord,Err)
    1529. 

  1529. ! ! !       CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,NODE,3,ZCoord,Err)
    1530. 

  1530. ! ! !       WRITE(*,*) "XCoord, YCoord, ZCoord = ",XCoord, YCoord, ZCoord
    1531. 

  1531. ! ! !     ENDIF
    1532. 

  1532. !     ENDDO
    1533. 

  1533. 

  1534.     !All other faces are impermeable
    1535. 

  1535.     DO NN=1,SIZE(Face7Nodes,1)
    1536. 

  1536.       NODE=Face7Nodes(NN)
    1537. 

  1537. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1538. 

  1538. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1539. 

  1539.         VALUE = 0.0_CMISSDP
    1540. 

  1540.         COMPONENT_NUMBER = 1
    1541. 

  1541.         write(*,*)'Marker 1'
    1542. 

  1542.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,1,1,NODE, &
    1543. 

  1543.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1544. 

  1544.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN X DIRECTION"
    1545. 

  1545. ! !     ENDIF
    1546. 

  1546.     ENDDO
    1547. 

  1547. 

  1548.     DO NN=1,SIZE(Face8Nodes,1)
    1549. 

  1549.       NODE=Face8Nodes(NN)
    1550. 

  1550. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1551. 

  1551. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1552. 

  1552.         VALUE = 0.0_CMISSDP
    1553. 

  1553.         COMPONENT_NUMBER = 1
    1554. 

  1554.         write(*,*)'Marker 2'
    1555. 

  1555.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,1,1,NODE, &
    1556. 

  1556.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1557. 

  1557.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN X DIRECTION"
    1558. 

  1558. ! !     ENDIF
    1559. 

  1559.     ENDDO
    1560. 

  1560. 

  1561.     DO NN=1,SIZE(Face9Nodes,1)
    1562. 

  1562.       NODE=Face9Nodes(NN)
    1563. 

  1563. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1564. 

  1564. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1565. 

  1565.         VALUE = 0.0_CMISSDP
    1566. 

  1566.         COMPONENT_NUMBER = 2
    1567. 

  1567.         write(*,*)'Marker 3'
    1568. 

  1568.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,1,1,NODE, &
    1569. 

  1569.           & COMPONENT_NUMBER, CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1570. 

  1570.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN Y DIRECTION"
    1571. 

  1571. ! !     ENDIF
    1572. 

  1572.     ENDDO
    1573. 

  1573. 

  1574.     DO NN=1,SIZE(Face10Nodes,1)
    1575. 

  1575.       NODE=Face10Nodes(NN)
    1576. 

  1576. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1577. 

  1577. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1578. 

  1578.         VALUE = 0.0_CMISSDP
    1579. 

  1579.         COMPONENT_NUMBER = 2
    1580. 

  1580.         write(*,*)'Marker 4'
    1581. 

  1581.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,1,1,NODE, &
    1582. 

  1582.           & COMPONENT_NUMBER,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1583. 

  1583.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN Y DIRECTION"
    1584. 

  1584. ! !     ENDIF
    1585. 

  1585.     ENDDO
    1586. 

  1586. 

  1587.     DO NN=1,SIZE(Face12Nodes,1)
    1588. 

  1588.       NODE=Face12Nodes(NN)
    1589. 

  1589. ! !     CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
    1590. 

  1590. ! !     IF(NodeDomain==ComputationalNodeNumber) THEN
    1591. 

  1591.         VALUE = 0.0_CMISSDP
    1592. 

  1592.         COMPONENT_NUMBER = 3
    1593. 

  1593.         write(*,*)'Marker 5'
    1594. 

  1594.         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,1,1,NODE, &
    1595. 

  1595.           & COMPONENT_NUMBER,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
    1596. 

  1596.         WRITE(*,*) "SPECIFIED IMPERMEABLE WALL AT NODE",NODE,"IN Z DIRECTION"
    1597. 

  1597. ! !     ENDIF
    1598. 

  1598.     ENDDO
    1599. 

  1599.     ENDIF
    1600. 

  1600. 

  1601.   !BOUNDARY CONDITIONS
    1602. 

  1602.   !Start the creation of the equations set boundary conditions for Darcy
    1603. 

  1603. !   DO icompartment=1,Ncompartments  
    1604. 

  1604. !     CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsDarcy,DependentFieldSolid,Err)
    1605. 

  1605. !     CALL CMISSEquationsSetBoundaryConditionsCreateStart(EquationsSetDarcy(icompartment),BoundaryConditionsDarcy,DependentFieldSolid,Err)
    1606. 

  1606. !     CALL CMISSEquationsSetBoundaryConditionsCreateFinish(EquationsSetDarcy(icompartment),Err)
    1607. 

  1607.   ENDDO
    1608. 

  1608.   CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsDarcy,Err)
    1609. 

  1609. 

  1610.   !
    1611. 

  1611.   !================================================================================================================================
    1612. 

  1612.   !
    1613. 

  1613. 

  1614.   !RUN SOLVERS
    1615. 

  1615. 

  1616.   !Turn of PETSc error handling
    1617. 

  1617.   !CALL PETSC_ERRORHANDLING_SET_ON(ERR,ERROR,*999)
    1618. 

  1618. 

  1619.   !Solve the problem
    1620. 

  1620.   WRITE(*,'(A)') "Solving problem..."
    1621. 

  1621.   CALL CMISSProblem_Solve(Problem,Err)
    1622. 

  1622.   WRITE(*,'(A)') "Problem solved!"
    1623. 

  1623. 

  1624. 

  1625.   !
    1626. 

  1626.   !================================================================================================================================
    1627. 

  1627.   !
    1628. 

  1628. 

  1629.   !OUTPUT
    1630. 

  1630. 

  1631.   WRITE(*,*)'dummy_counter = ',dummy_counter
    1632. 

  1632. 

  1633.   EXPORT_FIELD_IO=.FALSE.
    1634. 

  1634.   IF(EXPORT_FIELD_IO) THEN
    1635. 

  1635.     WRITE(*,'(A)') "Exporting fields..."
    1636. 

  1636.     CALL CMISSFields_Initialise(Fields,Err)
    1637. 

  1637.     CALL CMISSFields_Create(Region,Fields,Err)
    1638. 

  1638.     CALL CMISSFields_NodesExport(Fields,"FiniteElasticityMultiCompDarcy","FORTRAN",Err)
    1639. 

  1639.     CALL CMISSFields_ElementsExport(Fields,"FiniteElasticityMultiCompDarcy","FORTRAN",Err)
    1640. 

  1640.     CALL CMISSFields_Finalise(Fields,Err)
    1641. 

  1641.     WRITE(*,'(A)') "Field exported!"
    1642. 

  1642.   ENDIF
    1643. 

  1643. 

  1644. 

  1645.   !Finialise CMISS
    1646. 

  1646. !   CALL CMISSFinalise(Err)
    1647. 

  1647. 

  1648.   IF (ALLOCATED(EquationsSetFieldDarcy)) DEALLOCATE(EquationsSetFieldDarcy)
    1649. 

  1649.   IF (ALLOCATED(EquationsSetDarcy)) DEALLOCATE(EquationsSetDarcy)
    1650. 

  1650.   IF (ALLOCATED(MaterialsFieldDarcy)) DEALLOCATE(MaterialsFieldDarcy)
    1651. 

  1651.   IF (ALLOCATED(BoundaryConditionsDarcy)) DEALLOCATE(BoundaryConditionsDarcy)
    1652. 

  1652.   IF (ALLOCATED(EquationsDarcy)) DEALLOCATE(EquationsDarcy)
    1653. 

  1653. 

  1654.   WRITE(*,'(A)') "Program successfully completed."
    1655. 

  1655. 

  1656.   STOP
    1657. 

  1657. 

  1658. END PROGRAM FINITEELASTICITYMULTICOMPDARCYEXAMPLE
    1659. 

  1659.