PageRenderTime 140ms CodeModel.GetById 32ms app.highlight 97ms RepoModel.GetById 1ms app.codeStats 1ms

/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/IncompElastDrivenDarcy_AnalyticDarcy/src/IncompElastDrivenDarcy_AnalyticDarcyExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 1379 lines | 566 code | 185 blank | 628 comment | 0 complexity | 11cc114a73868b9bf78a2c085958a0ff MD5 | raw file
   1!> \file
   2!> \author Christian Michler, Adam Reeve
   3!> \brief This is an example program to solve a coupled Finite Elastiticity Darcy equation using openCMISS calls.
   4!>
   5!> \section LICENSE
   6!>
   7!> Version: MPL 1.1/GPL 2.0/LGPL 2.1
   8!>
   9!> The contents of this file are subject to the Mozilla Public License
  10!> Version 1.1 (the "License"); you may not use this file except in
  11!> compliance with the License. You may obtain a copy of the License at
  12!> http://www.mozilla.org/MPL/
  13!>
  14!> Software distributed under the License is distributed on an "AS IS"
  15!> basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
  16!> License for the specific language governing rights and limitations
  17!> under the License.
  18!>
  19!> The Original Code is OpenCMISS
  20!>
  21!> The Initial Developer of the Original Code is University of Auckland,
  22!> Auckland, New Zealand and University of Oxford, Oxford, United
  23!> Kingdom. Portions created by the University of Auckland and University
  24!> of Oxford are Copyright (C) 2007 by the University of Auckland and
  25!> the University of Oxford. All Rights Reserved.
  26!>
  27!> Contributor(s):
  28!>
  29!> Alternatively, the contents of this file may be used under the terms of
  30!> either the GNU General Public License Version 2 or later (the "GPL"), or
  31!> the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
  32!> in which case the provisions of the GPL or the LGPL are applicable instead
  33!> of those above. If you wish to allow use of your version of this file only
  34!> under the terms of either the GPL or the LGPL, and not to allow others to
  35!> use your version of this file under the terms of the MPL, indicate your
  36!> decision by deleting the provisions above and replace them with the notice
  37!> and other provisions required by the GPL or the LGPL. If you do not delete
  38!> the provisions above, a recipient may use your version of this file under
  39!> the terms of any one of the MPL, the GPL or the LGPL.
  40!>
  41
  42!> \example MultiPhysics/Poroelasticity/FiniteElasticityDarcy/IncompressibleElasticityDrivenDarcy/src/IncompressibleElasticityDrivenDarcyExample.f90
  43!! Example program to solve coupled FiniteElasticityDarcy equations using OpenCMISS calls.
  44!! \par Latest Builds:
  45!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/IncompressibleElasticityDrivenDarcy/build-intel'>Linux Intel Build</a>
  46!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/IncompressibleElasticityDrivenDarcy/build-intel'>Linux GNU Build</a>
  47!!
  48!<
  49
  50! !
  51! !  This example considers a coupled Finite Elasticity Darcy problem
  52! !
  53
  54!> Main program
  55
  56PROGRAM INCOMPELASTDRIVENDARCYANALYTICDARCYEXAMPLE
  57
  58  !
  59  !================================================================================================================================
  60  !
  61
  62  !PROGRAM LIBRARIES
  63
  64  USE OPENCMISS
  65  USE FLUID_MECHANICS_IO_ROUTINES
  66  USE MPI
  67
  68#ifdef WIN32
  69  USE IFQWINCMISS
  70#endif
  71
  72  !
  73  !================================================================================================================================
  74  !
  75
  76  !PROGRAM VARIABLES AND TYPES
  77
  78  IMPLICIT NONE
  79
  80  !Test program parameters
  81
  82  REAL(CMISSDP), PARAMETER :: HEIGHT=1.0_CMISSDP
  83  REAL(CMISSDP), PARAMETER :: WIDTH=1.0_CMISSDP
  84  REAL(CMISSDP), PARAMETER :: LENGTH=1.0_CMISSDP
  85
  86  INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=1
  87  INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=2
  88  INTEGER(CMISSIntg), PARAMETER :: CubicBasisUserNumber=3
  89
  90  INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
  91  INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=2
  92  INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=3
  93  INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=4
  94  INTEGER(CMISSIntg), PARAMETER :: GeometricFieldUserNumber=5
  95  INTEGER(CMISSIntg), PARAMETER :: DependentFieldUserNumberDarcy=6
  96  INTEGER(CMISSIntg), PARAMETER :: DependentFieldUserNumberMatProperties=42
  97  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcy=8
  98  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberMatProperties=9
  99  INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberDarcy=12
 100  INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberMatProperties=13
 101  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=14
 102  INTEGER(CMISSIntg), PARAMETER :: IndependentFieldUserNumberSolid=15
 103  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumberDarcy=22
 104  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumberMatProperties=23
 105  INTEGER(CMISSIntg), PARAMETER :: AnalyticFieldUserNumberDarcy=27
 106
 107
 108  INTEGER(CMISSIntg), PARAMETER :: NumberOfDomains=1
 109  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1
 110  INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2
 111  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1
 112  INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1
 113  INTEGER(CMISSIntg), PARAMETER :: SolverMatPropertiesIndex=1
 114  INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=2
 115  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPorosity=1
 116  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPermOverVis=2
 117  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberMatPropertiesPorosity=1
 118  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberMatPropertiesPermOverVis=2
 119
 120  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1
 121  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
 122
 123  !Program types
 124
 125  TYPE(EXPORT_CONTAINER):: CM
 126
 127  !Program variables
 128
 129  INTEGER(CMISSIntg) :: NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS,NUMBER_GLOBAL_Z_ELEMENTS
 130  INTEGER(CMISSIntg) :: NUMBER_OF_DOMAINS
 131
 132  INTEGER(CMISSIntg) :: MPI_IERROR
 133
 134  INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
 135
 136  INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
 137  INTEGER(CMISSIntg) :: RESTART_VALUE
 138
 139  INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT
 140  INTEGER(CMISSIntg) :: EQUATIONS_MAT_PROPERTIES_OUTPUT
 141  INTEGER(CMISSIntg) :: COMPONENT_NUMBER
 142  INTEGER(CMISSIntg) :: NODE_NUMBER
 143  INTEGER(CMISSIntg) :: ELEMENT_NUMBER
 144  INTEGER(CMISSIntg) :: CONDITION
 145
 146  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY
 147  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE
 148  INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE
 149  INTEGER(CMISSIntg) :: LINEAR_SOLVER_MAT_PROPERTIES_OUTPUT_TYPE
 150
 151  REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z
 152  REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2
 153  REAL(CMISSDP) :: GEOMETRY_TOLERANCE
 154  INTEGER(CMISSIntg) :: EDGE_COUNT
 155  INTEGER(CMISSIntg) :: NUMBER_OF_COMPONENTS_DEPENDENT_FIELD_MAT_PROPERTIES
 156  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID
 157  REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4)
 158  REAL(CMISSDP) :: INITIAL_FIELD_MAT_PROPERTIES(3)
 159  REAL(CMISSDP) :: INITIAL_FIELD_SOLID(4)
 160  REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
 161  REAL(CMISSDP) :: RELATIVE_TOLERANCE
 162  REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
 163  REAL(CMISSDP) :: LINESEARCH_ALPHA
 164  REAL(CMISSDP) :: VALUE
 165  REAL(CMISSDP) :: POROSITY_PARAM_MAT_PROPERTIES, PERM_OVER_VIS_PARAM_MAT_PROPERTIES
 166  REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY
 167
 168  LOGICAL :: EXPORT_FIELD_IO
 169  LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG
 170  LOGICAL :: LINEAR_SOLVER_MAT_PROPERTIES_DIRECT_FLAG
 171
 172  !CMISS variables
 173
 174  !Regions
 175  TYPE(CMISSRegionType) :: Region
 176  TYPE(CMISSRegionType) :: WorldRegion
 177  !Coordinate systems
 178  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem
 179  TYPE(CMISSCoordinateSystemType) :: WorldCoordinateSystem
 180  !Basis
 181  TYPE(CMISSBasisType) :: BasisGeometry
 182  TYPE(CMISSBasisType) :: BasisVelocity
 183  TYPE(CMISSBasisType) :: BasisPressure
 184  TYPE(CMISSBasisType) :: CubicBasis, QuadraticBasis, LinearBasis, Bases(2)
 185  !Nodes
 186  TYPE(CMISSNodesType) :: Nodes
 187  !Elements
 188  TYPE(CMISSMeshElementsType) :: MeshElementsGeometry
 189  TYPE(CMISSMeshElementsType) :: MeshElementsVelocity
 190  TYPE(CMISSMeshElementsType) :: MeshElementsPressure
 191  !Meshes
 192  TYPE(CMISSMeshType) :: Mesh
 193  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
 194
 195  !Decompositions
 196  TYPE(CMISSDecompositionType) :: Decomposition
 197  !Fields
 198  TYPE(CMISSFieldsType) :: Fields
 199  !Field types
 200  TYPE(CMISSFieldType) :: GeometricField
 201  TYPE(CMISSFieldType) :: DependentFieldMatProperties
 202  TYPE(CMISSFieldType) :: MaterialsFieldDarcy
 203  TYPE(CMISSFieldType) :: MaterialsFieldMatProperties
 204  TYPE(CMISSFieldType) :: EquationsSetFieldDarcy
 205  TYPE(CMISSFieldType) :: EquationsSetFieldMatProperties
 206  TYPE(CMISSFieldType) :: IndependentFieldSolid
 207  TYPE(CMISSFieldType) :: AnalyticFieldDarcy
 208  !Boundary conditions
 209  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy
 210  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsMatProperties
 211  !Equations sets
 212  TYPE(CMISSEquationsSetType) :: EquationsSetDarcy
 213  TYPE(CMISSEquationsSetType) :: EquationsSetMatProperties
 214  !Equations
 215  TYPE(CMISSEquationsType) :: EquationsDarcy
 216  TYPE(CMISSEquationsType) :: EquationsMatProperties
 217  !Problems
 218  TYPE(CMISSProblemType) :: Problem
 219  !Control loops
 220  TYPE(CMISSControlLoopType) :: ControlLoop
 221  !Solvers
 222  TYPE(CMISSSolverType) :: DynamicSolverDarcy
 223  TYPE(CMISSSolverType) :: LinearSolverDarcy
 224  TYPE(CMISSSolverType) :: LinearSolverMatProperties
 225  TYPE(CMISSSolverType) :: LinearSolverSolid
 226  !Solver equations
 227  TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy
 228  TYPE(CMISSSolverEquationsType) :: SolverEquationsMatProperties
 229
 230#ifdef WIN32
 231  !Quickwin type
 232  LOGICAL :: QUICKWIN_STATUS=.FALSE.
 233  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
 234#endif
 235
 236  !Generic CMISS variables
 237
 238  INTEGER(CMISSIntg) :: EquationsSetIndex
 239  INTEGER(CMISSIntg) :: Err
 240
 241
 242  INTEGER(CMISSIntg) :: DIAG_LEVEL_LIST(5)
 243!   CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(8) !,TIMING_ROUTINE_LIST(1)
 244  CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(1) !,TIMING_ROUTINE_LIST(1)
 245
 246
 247
 248  !
 249  !--------------------------------------------------------------------------------------------------------------------------------
 250  !
 251
 252  !Program variables and types (finite elasticity part)
 253
 254  !Test program parameters
 255
 256  INTEGER(CMISSIntg) :: BASIS_NUMBER_SOLID
 257
 258  INTEGER(CMISSIntg) :: TotalNumberOfSolidNodes
 259  INTEGER(CMISSIntg) :: SolidMeshComponenetNumber
 260  INTEGER(CMISSIntg) :: SolidPressureMeshComponenetNumber
 261
 262  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidUserNumber=1
 263  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfVariables=1
 264  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfComponents=3
 265
 266  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidUserNumber=2
 267  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfVariables=1
 268  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfComponents=3
 269
 270  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidUserNumber=3
 271  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfVariables=1
 272  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfComponents=3
 273
 274  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidUserNumber=4
 275  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfVariables=4
 276  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4
 277  INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4  !(u,v,w,m)
 278
 279  INTEGER(CMISSIntg), PARAMETER :: EquationSetSolidUserNumber=1
 280  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldSolidUserNumber=25
 281
 282  INTEGER(CMISSIntg), PARAMETER :: DisplacementMeshComponentNumber=1
 283  INTEGER(CMISSIntg), PARAMETER :: PressureMeshComponentNumber=2
 284
 285  INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=32
 286  !Program types
 287  !Program variables
 288
 289  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME
 290  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME
 291  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA
 292  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
 293
 294  !CMISS variables
 295
 296  TYPE(CMISSBasisType) :: BasisSolid
 297  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsSolid
 298  TYPE(CMISSEquationsType) :: EquationsSolid
 299  TYPE(CMISSEquationsSetType) :: EquationsSetSolid
 300  TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid
 301  TYPE(CMISSFieldType) :: DependentFieldSolid,EquationsSetFieldSolid
 302  TYPE(CMISSSolverType) :: SolverSolid
 303  TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid
 304  TYPE(CMISSMeshElementsType) :: MeshElementsSolid
 305
 306  !End - Program variables and types (finite elasticity part)
 307
 308  !
 309  !--------------------------------------------------------------------------------------------------------------------------------
 310  !
 311
 312
 313#ifdef WIN32
 314  !Initialise QuickWin
 315  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
 316  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
 317  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
 318  !Set the window parameters
 319  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
 320  !If attempt fails set with system estimated values
 321  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
 322#endif
 323
 324  !
 325  !================================================================================================================================
 326  !
 327  NUMBER_GLOBAL_X_ELEMENTS=2
 328  NUMBER_GLOBAL_Y_ELEMENTS=2
 329  NUMBER_GLOBAL_Z_ELEMENTS=2
 330
 331  IF(NUMBER_GLOBAL_Z_ELEMENTS==0)THEN
 332    NUMBER_OF_DIMENSIONS=2
 333  ELSE
 334    NUMBER_OF_DIMENSIONS=3
 335  ENDIF
 336  !PROBLEM CONTROL PANEL
 337
 338  !Import cmHeart mesh information
 339  !CALL FLUID_MECHANICS_IO_READ_CMHEART(CM,Err)
 340  BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION
 341  !Set geometric tolerance
 342  GEOMETRY_TOLERANCE = 1.0E-12_CMISSDP
 343  !Set initial values
 344  INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
 345  INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
 346  INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
 347  INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
 348  INITIAL_FIELD_MAT_PROPERTIES(1)=0.0_CMISSDP
 349  INITIAL_FIELD_MAT_PROPERTIES(2)=0.0_CMISSDP
 350  INITIAL_FIELD_MAT_PROPERTIES(3)=0.0_CMISSDP
 351!   INITIAL_FIELD_SOLID(1)=1.0_CMISSDP
 352!   INITIAL_FIELD_SOLID(2)=1.0_CMISSDP
 353!   INITIAL_FIELD_SOLID(3)=1.0_CMISSDP
 354!   INITIAL_FIELD_SOLID(4)=1.0_CMISSDP
 355  !Set material parameters
 356  POROSITY_PARAM_DARCY=0.1_CMISSDP
 357  PERM_OVER_VIS_PARAM_DARCY=1.0e-1_CMISSDP
 358!   PERM_OVER_VIS_PARAM_DARCY=0.1_CMISSDP
 359  POROSITY_PARAM_MAT_PROPERTIES=POROSITY_PARAM_DARCY
 360  PERM_OVER_VIS_PARAM_MAT_PROPERTIES=PERM_OVER_VIS_PARAM_DARCY
 361  !Set output parameter
 362  !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
 363  LINEAR_SOLVER_MAT_PROPERTIES_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
 364  DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
 365  LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_SOLVER_OUTPUT
 366  !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
 367  EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
 368  EQUATIONS_MAT_PROPERTIES_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
 369
 370  !Set time parameter
 371  DYNAMIC_SOLVER_DARCY_START_TIME=0.0_CMISSDP
 372!   DYNAMIC_SOLVER_DARCY_STOP_TIME=0.03_CMISSDP
 373  DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-2_CMISSDP
 374  DYNAMIC_SOLVER_DARCY_STOP_TIME=2_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
 375  DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
 376  !Set result output parameter
 377  DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
 378  !Set solver parameters
 379  LINEAR_SOLVER_MAT_PROPERTIES_DIRECT_FLAG=.TRUE.
 380  LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
 381  RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
 382  ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
 383  DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
 384  MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
 385  RESTART_VALUE=30_CMISSIntg !default: 30
 386  LINESEARCH_ALPHA=1.0_CMISSDP
 387
 388
 389  !
 390  !================================================================================================================================
 391  !
 392
 393  !INITIALISE OPENCMISS
 394
 395  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
 396
 397  !CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
 398
 399  !
 400  !================================================================================================================================
 401  !
 402
 403  !Set diagnostics
 404
 405  DIAG_LEVEL_LIST(1)=1
 406  DIAG_LEVEL_LIST(2)=2
 407  DIAG_LEVEL_LIST(3)=3
 408  DIAG_LEVEL_LIST(4)=4
 409  DIAG_LEVEL_LIST(5)=5
 410
 411!   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_FINITE_ELEMENT_CALCULATE"
 412!   DIAG_ROUTINE_LIST(2)="DARCY_EQUATION_PRE_SOLVE_STORE_REFERENCE_DATA"
 413!   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_GET_SOLID_DISPLACEMENT"
 414!   DIAG_ROUTINE_LIST(2)="DARCY_EQUATION_PRE_SOLVE_ALE_UPDATE_MESH"
 415!   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_GET_SOLID_DISPLACEMENT"
 416!   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_PRE_SOLVE_UPDATE_BOUNDARY_CONDITIONS"
 417!   DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
 418!   DIAG_ROUTINE_LIST(1)="DARCY_EQUATION_POST_SOLVE_ADD_MASS_CORRECTION"
 419  DIAG_ROUTINE_LIST(1)="WRITE_IP_INFO"
 420!   DIAG_ROUTINE_LIST(2)="FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR"
 421!   DIAG_ROUTINE_LIST(3)="EVALUATE_CHAPELLE_PIOLA_TENSOR_ADDITION"
 422!   DIAG_ROUTINE_LIST(5)="DARCY_EQUATION_PRE_SOLVE_MAT_PROPERTIES"
 423!   DIAG_ROUTINE_LIST(6)="FITTING_FINITE_ELEMENT_CALCULATE"
 424!   DIAG_ROUTINE_LIST(7)="FINITE_ELASTICITY_FINITE_ELEMENT_JACOBIAN_EVALUATE"
 425!   DIAG_ROUTINE_LIST(8)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
 426!   DIAG_ROUTINE_LIST(1)="PROBLEM_SOLVER_EQUATIONS_SOLVE"
 427!   DIAG_ROUTINE_LIST(1)="SOLVER_NEWTON_SOLVE"
 428!   DIAG_ROUTINE_LIST(2)="SOLVER_NEWTON_LINESEARCH_SOLVE"
 429!   DIAG_ROUTINE_LIST(1)="SOLVER_SOLUTION_UPDATE"
 430!   DIAG_ROUTINE_LIST(1)="FINITE_ELASTICITY_FINITE_ELEMENT_RESIDUAL_EVALUATE"
 431
 432  !CMISS_ALL_DIAG_TYPE/CMISS_IN_DIAG_TYPE/CMISS_FROM_DIAG_TYPE
 433  CALL CMISSDiagnosticsSetOn(CMISS_IN_DIAG_TYPE,DIAG_LEVEL_LIST,"Diagnostics",DIAG_ROUTINE_LIST,Err)
 434
 435  !CMISS_ALL_TIMING_TYPE/CMISS_IN_TIMING_TYPE/CMISS_FROM_TIMING_TYPE
 436  !TIMING_ROUTINE_LIST(1)="PROBLEM_FINITE_ELEMENT_CALCULATE"
 437  !CALL TIMING_SET_ON(IN_TIMING_TYPE,.TRUE.,"",TIMING_ROUTINE_LIST,ERR,ERROR,*999)
 438
 439  !
 440  !================================================================================================================================
 441  !
 442
 443  !COORDINATE SYSTEM
 444
 445  !Start the creation of a new RC coordinate system
 446  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
 447  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
 448  !Set the coordinate system dimension
 449  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NUMBER_OF_DIMENSIONS,Err)
 450  !Finish the creation of the coordinate system
 451  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
 452
 453  !
 454  !================================================================================================================================
 455  !
 456
 457  !REGION
 458  !For a volume-coupled problem, solid and fluid are based in the same region
 459
 460  !Start the creation of a new region
 461  CALL CMISSRegion_Initialise(Region,Err)
 462  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
 463  !Set the regions coordinate system as defined above
 464  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
 465  !Finish the creation of the region
 466  CALL CMISSRegion_CreateFinish(Region,Err)
 467
 468  !
 469  !================================================================================================================================
 470  !
 471
 472  !BASES
 473  !Define basis functions
 474  CALL CMISSBasis_Initialise(LinearBasis,Err)
 475  CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
 476  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
 477    & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
 478  !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err)
 479  CALL CMISSBasis_CreateFinish(LinearBasis,Err)
 480
 481  CALL CMISSBasis_Initialise(QuadraticBasis,Err)
 482  CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
 483  CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
 484    & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
 485  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
 486    & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
 487  !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err)
 488  CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
 489
 490  CALL CMISSBasis_Initialise(CubicBasis,Err)
 491  CALL CMISSBasis_CreateStart(CubicBasisUserNumber,CubicBasis,Err)
 492  CALL CMISSBasis_InterpolationXiSet(CubicBasis,(/CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION, &
 493    & CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION/),Err)
 494  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(CubicBasis, &
 495    & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
 496  !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(CubicBasis,.true.,Err) !Enable 3D interpolation on faces
 497  CALL CMISSBasis_CreateFinish(CubicBasis,Err)
 498
 499  Bases(1)=CubicBasis
 500  Bases(2)=QuadraticBasis
 501
 502  !Start the creation of a generated mesh in the region
 503  CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
 504  CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
 505  !Set up a regular x*y*z mesh
 506  CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_REGULAR_MESH_TYPE,Err)
 507  !Set the default basis
 508  !CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,BasisGeometry,Err)
 509  CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,Bases,Err)
 510  !Define the mesh on the region
 511  IF(NUMBER_GLOBAL_Z_ELEMENTS==0) THEN
 512    CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/WIDTH,HEIGHT/),Err)
 513    CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS/),Err)
 514  ELSE
 515    CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/WIDTH,HEIGHT,LENGTH/),Err)
 516    CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS, &
 517      & NUMBER_GLOBAL_Z_ELEMENTS/),Err)
 518  ENDIF
 519  !Finish the creation of a generated mesh in the region
 520  CALL CMISSMesh_Initialise(Mesh,Err)
 521  CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
 522
 523  !GEOMETRIC FIELD
 524
 525  !Create a decomposition:
 526  !All mesh components (associated with G.Projection / Darcy / solid) share the same decomposition
 527  CALL CMISSDecomposition_Initialise(Decomposition,Err)
 528  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
 529  !Set the decomposition to be a general decomposition with the specified number of domains
 530  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
 531  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
 532  !Finish the decomposition
 533  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
 534
 535  !Start to create a default (geometric) field on the region
 536  CALL CMISSField_Initialise(GeometricField,Err)
 537  CALL CMISSField_CreateStart(GeometricFieldUserNumber,Region,GeometricField,Err)
 538  CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
 539  !Set the field type
 540  CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
 541  !Set the decomposition to use
 542  CALL CMISSField_NumberOfVariablesSet(GeometricField,FieldGeometryNumberOfVariables,Err)
 543  CALL CMISSField_NumberOfComponentsSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,Err)  
 544  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,1,DisplacementMeshComponentNumber,Err)
 545  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,2,DisplacementMeshComponentNumber,Err)
 546  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,DisplacementMeshComponentNumber,Err)
 547  CALL CMISSField_CreateFinish(GeometricField,Err)
 548  !Set the mesh component to be used by the field components.
 549  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricField,Err)
 550
 551!   !Update the geometric field parameters
 552!   DO NODE_NUMBER=1,NUMBER_OF_NODES_GEOMETRY
 553!     DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 554!       VALUE=CM%N(NODE_NUMBER,COMPONENT_NUMBER)
 555!       CALL CMISSField_ParameterSetUpdateNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 556!         & CMISS_NO_GLOBAL_DERIV,NODE_NUMBER,COMPONENT_NUMBER,VALUE,Err)
 557!     ENDDO
 558!   ENDDO
 559!   CALL CMISSField_ParameterSetUpdateStart(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
 560!   CALL CMISSField_ParameterSetUpdateFinish(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
 561
 562  !--------------------------------------------------------------------------------------------------------------------------------
 563  ! Solid
 564
 565  !Create a decomposition
 566
 567!   !Create a field to put the geometry (defualt is geometry)
 568!   CALL CMISSField_Initialise(GeometricFieldSolid,Err)
 569!   CALL CMISSField_CreateStart(FieldGeometrySolidUserNumber,Region,GeometricFieldSolid,Err)
 570!   CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err)
 571!   CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err)
 572!   CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometrySolidNumberOfVariables,Err)
 573!   CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometrySolidNumberOfComponents,Err)
 574!   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidMeshComponenetNumber,Err)
 575!   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidMeshComponenetNumber,Err)
 576!   CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidMeshComponenetNumber,Err)
 577!   CALL CMISSField_CreateFinish(GeometricFieldSolid,Err)
 578! 
 579! !---
 580!   !Update the geometric field parameters
 581!   DO NODE_NUMBER=1,NUMBER_OF_NODES_GEOMETRY
 582!     DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 583!       VALUE=CM%N(NODE_NUMBER,COMPONENT_NUMBER)
 584!       CALL CMISSField_ParameterSetUpdateNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 585!         & CMISS_NO_GLOBAL_DERIV,NODE_NUMBER,COMPONENT_NUMBER,VALUE,Err)
 586!     ENDDO
 587!   ENDDO
 588!   CALL CMISSField_ParameterSetUpdateStart(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
 589!   CALL CMISSField_ParameterSetUpdateFinish(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
 590!---
 591
 592  !Create a fibre field and attach it to the geometric field
 593  CALL CMISSField_Initialise(FibreFieldSolid,Err)
 594  CALL CMISSField_CreateStart(FieldFibreSolidUserNumber,Region,FibreFieldSolid,Err)
 595  CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err)
 596  CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err)
 597  CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricField,Err)
 598  CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreSolidNumberOfVariables,Err)
 599  CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreSolidNumberOfComponents,Err)
 600  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,DisplacementMeshComponentNumber,Err)
 601  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,DisplacementMeshComponentNumber,Err)
 602  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,DisplacementMeshComponentNumber,Err)
 603  CALL CMISSField_CreateFinish(FibreFieldSolid,Err)
 604
 605  ! end Solid
 606  !--------------------------------------------------------------------------------------------------------------------------------
 607
 608  !
 609  !================================================================================================================================
 610  !
 611
 612  !EQUATIONS SETS
 613
 614  !Create the equations set for ALE Darcy
 615  CALL CMISSField_Initialise(EquationsSetFieldDarcy,Err)
 616  CALL CMISSEquationsSet_Initialise(EquationsSetDarcy,Err)
 617  CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricField,CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
 618    & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,&
 619    & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy,EquationsSetDarcy,Err)
 620  !Set the equations set to be a ALE Darcy problem
 621!   CALL CMISSEquationsSet_SpecificationSet(EquationsSetDarcy,CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
 622!     & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,Err)
 623  !Finish creating the equations set
 624  CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy,Err)
 625
 626  !Create the equations set for deformation-dependent material properties
 627  CALL CMISSField_Initialise(EquationsSetFieldMatProperties,Err)
 628  CALL CMISSEquationsSet_Initialise(EquationsSetMatProperties,Err)
 629!   CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberMatProperties,Region,GeometricField,CMISS_EQUATIONS_SET_FITTING_CLASS,&
 630  CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberMatProperties,Region,GeometricField,CMISS_EQUATIONS_SET_FITTING_CLASS,&
 631    & CMISS_EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE,CMISS_EQUATIONS_SET_MAT_PROP_INRIA_MODEL_DATA_FITTING_SUBTYPE,&
 632    & EquationsSetFieldUserNumberMatProperties,EquationsSetFieldMatProperties,EquationsSetMatProperties,Err)
 633  !Set the equations set to be a deformation-dependent material properties problem
 634!   CALL CMISSEquationsSet_SpecificationSet(EquationsSetMatProperties,CMISS_EQUATIONS_SET_FITTING_CLASS, &
 635!     & CMISS_EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE,CMISS_EQUATIONS_SET_MAT_PROP_INRIA_MODEL_DATA_FITTING_SUBTYPE,Err)
 636  !Finish creating the equations set
 637  CALL CMISSEquationsSet_CreateFinish(EquationsSetMatProperties,Err)
 638
 639
 640  !--------------------------------------------------------------------------------------------------------------------------------
 641  ! Solid
 642
 643  !Create the equations_set
 644  CALL CMISSField_Initialise(EquationsSetFieldSolid,Err)
 645  CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err)
 646  CALL CMISSEquationsSet_CreateStart(EquationSetSolidUserNumber,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
 647    & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,&
 648    & EquationsSetFieldSolidUserNumber,EquationsSetFieldSolid,EquationsSetSolid,Err)
 649!   CALL CMISSEquationsSet_SpecificationSet(EquationsSetSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
 650!     & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,Err)
 651  CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err)
 652
 653  ! end Solid
 654  !--------------------------------------------------------------------------------------------------------------------------------
 655  !--------------------------------------------------------------------------------------------------------------------------------
 656  ! Solid Materials Field
 657
 658  !Create a material field and attach it to the geometric field
 659  CALL CMISSField_Initialise(MaterialFieldSolid,Err)
 660  !
 661  CALL CMISSField_CreateStart(FieldMaterialSolidUserNumber,Region,MaterialFieldSolid,Err)
 662  !
 663  CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err)
 664  CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err)
 665  CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricField,Err)
 666  CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialSolidNumberOfVariables,Err)
 667  CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialSolidNumberOfComponents,Err)
 668  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,DisplacementMeshComponentNumber,Err)
 669  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,DisplacementMeshComponentNumber,Err)
 670  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,DisplacementMeshComponentNumber,Err)
 671  !
 672  CALL CMISSField_CreateFinish(MaterialFieldSolid,Err)
 673
 674  !Set material parameters
 675  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
 676    & 1.0_CMISSDP,Err)
 677!   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,2.0e3_CMISSDP,Err)
 678  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2, &
 679    & 1.0_CMISSDP,Err)
 680!   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,33.0_CMISSDP,Err)
 681  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3, &
 682    & 10.0_CMISSDP,Err)
 683
 684  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialSolidUserNumber,MaterialFieldSolid,Err)
 685  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err)
 686
 687  ! end Solid
 688  !--------------------------------------------------------------------------------------------------------------------------------
 689
 690
 691  !
 692  !================================================================================================================================
 693  !
 694
 695  !DEPENDENT FIELDS
 696
 697  !Create the equations set dependent field variables for deformation-dependent material properties
 698  CALL CMISSField_Initialise(DependentFieldMatProperties,Err)
 699  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetMatProperties,DependentFieldUserNumberMatProperties, &
 700    & DependentFieldMatProperties,Err)
 701  !Set the mesh component to be used by the field components.
 702  NUMBER_OF_COMPONENTS_DEPENDENT_FIELD_MAT_PROPERTIES = 2
 703  DO COMPONENT_NUMBER=1,NUMBER_OF_COMPONENTS_DEPENDENT_FIELD_MAT_PROPERTIES
 704    CALL CMISSField_ComponentMeshComponentSet(DependentFieldMatProperties,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
 705      & DisplacementMeshComponentNumber,Err)
 706    CALL CMISSField_ComponentMeshComponentSet(DependentFieldMatProperties,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
 707      & DisplacementMeshComponentNumber,Err)
 708  ENDDO
 709  !Finish the equations set dependent field variables
 710  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetMatProperties,Err)
 711
 712  !Initialise dependent field
 713  DO COMPONENT_NUMBER=1,NUMBER_OF_COMPONENTS_DEPENDENT_FIELD_MAT_PROPERTIES
 714    CALL CMISSField_ComponentValuesInitialise(DependentFieldMatProperties,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 715      & COMPONENT_NUMBER,INITIAL_FIELD_MAT_PROPERTIES(COMPONENT_NUMBER),Err)
 716  ENDDO
 717
 718  !--------------------------------------------------------------------------------------------------------------------------------
 719  ! Solid
 720
 721  !Create a dependent field with two variables and four components
 722  CALL CMISSField_Initialise(DependentFieldSolid,Err)
 723  !
 724  CALL CMISSField_CreateStart(FieldDependentSolidUserNumber,Region,DependentFieldSolid,Err)
 725  !
 726  CALL CMISSField_TypeSet(DependentFieldSolid,CMISS_FIELD_GENERAL_TYPE,Err)
 727  CALL CMISSField_MeshDecompositionSet(DependentFieldSolid,Decomposition,Err)
 728  CALL CMISSField_GeometricFieldSet(DependentFieldSolid,GeometricField,Err)
 729  CALL CMISSField_DependentTypeSet(DependentFieldSolid,CMISS_FIELD_DEPENDENT_TYPE,Err)
 730  CALL CMISSField_NumberOfVariablesSet(DependentFieldSolid,FieldDependentSolidNumberOfVariables,Err)
 731  CALL CMISSField_VariableTypesSet(DependentFieldSolid,(/CMISS_FIELD_U_VARIABLE_TYPE, &
 732    & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err)
 733  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
 734  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
 735    & FieldDependentSolidNumberOfComponents,Err)
 736  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
 737  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE, &
 738    & FieldDependentFluidNumberOfComponents,Err)
 739  !
 740  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,DisplacementMeshComponentNumber,Err)
 741  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,DisplacementMeshComponentNumber,Err)
 742  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,DisplacementMeshComponentNumber,Err)
 743  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4, &
 744    & CMISS_FIELD_NODE_BASED_INTERPOLATION, &
 745    & Err)
 746!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
 747!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
 748  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,PressureMeshComponentNumber,Err)
 749  !
 750  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
 751    & DisplacementMeshComponentNumber,Err)
 752  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
 753    & DisplacementMeshComponentNumber,Err)
 754  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
 755    & DisplacementMeshComponentNumber,Err)
 756  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
 757    & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
 758!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
 759!     & CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
 760!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
 761  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,PressureMeshComponentNumber, &
 762    & Err)
 763
 764  !For this equation type, MESH_COMPONENT_NUMBER_PRESSURE is actually the mass increase component as the pressure is taken from the solid equations
 765  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,DisplacementMeshComponentNumber,Err)
 766  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,2,DisplacementMeshComponentNumber,Err)
 767  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,3,DisplacementMeshComponentNumber,Err)
 768!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
 769  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,PressureMeshComponentNumber,Err)
 770  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1, &
 771    & DisplacementMeshComponentNumber,Err)
 772  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2, &
 773    & DisplacementMeshComponentNumber,Err)
 774  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3, &
 775    & DisplacementMeshComponentNumber,Err)
 776!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
 777  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,PressureMeshComponentNumber, &
 778    & Err)
 779
 780  !
 781  CALL CMISSField_CreateFinish(DependentFieldSolid,Err)
 782  !
 783  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentSolidUserNumber,DependentFieldSolid,Err)
 784  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err)
 785
 786!   !Initialise dependent field (solid displacement and pressure)
 787!   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS !+1
 788!     CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 789!       & COMPONENT_NUMBER,INITIAL_FIELD_SOLID(COMPONENT_NUMBER),Err)
 790!   ENDDO
 791
 792  ! end Solid
 793  !--------------------------------------------------------------------------------------------------------------------------------
 794
 795
 796  !Create the equations set dependent field variables for ALE Darcy
 797!   CALL CMISSField_Initialise(DependentFieldDarcy,Err)
 798!   CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,DependentFieldUserNumberDarcy, & ! ??? UserNumber ???
 799  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,FieldDependentSolidUserNumber, & ! ??? UserNumber ???
 800    & DependentFieldSolid,Err)
 801!   !Set the mesh component to be used by the field components.
 802!   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 803!     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
 804!       & MESH_COMPONENT_NUMBER_VELOCITY,Err)
 805!     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
 806!       & MESH_COMPONENT_NUMBER_VELOCITY,Err)
 807!   ENDDO
 808!   COMPONENT_NUMBER=NUMBER_OF_DIMENSIONS+1
 809!     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
 810!       & MESH_COMPONENT_NUMBER_PRESSURE,Err)
 811!     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDarcy,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
 812!       & MESH_COMPONENT_NUMBER_PRESSURE,Err)
 813  !Finish the equations set dependent field variables
 814  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy,Err)
 815
 816  !Initialise dependent field (velocity components,pressure,mass increase)
 817  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS+1
 818    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 819      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
 820  ENDDO
 821
 822
 823  !
 824  !================================================================================================================================
 825  !
 826
 827  !MATERIALS FIELDS
 828
 829  !Create the equations set materials field variables for ALE Darcy
 830  CALL CMISSField_Initialise(MaterialsFieldDarcy,Err)
 831  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy,MaterialsFieldUserNumberDarcy, &
 832    & MaterialsFieldDarcy,Err)
 833  !Finish the equations set materials field variables
 834  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy,Err)
 835  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 836    & MaterialsFieldUserNumberDarcyPorosity,POROSITY_PARAM_DARCY,Err)
 837  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 838    & MaterialsFieldUserNumberDarcyPermOverVis,PERM_OVER_VIS_PARAM_DARCY,Err)
 839  !Create the equations set materials field variables for deformation-dependent material properties
 840  CALL CMISSField_Initialise(MaterialsFieldMatProperties,Err)
 841  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetMatProperties,MaterialsFieldUserNumberMatProperties, &
 842    & MaterialsFieldMatProperties,Err)
 843  !Finish the equations set materials field variables
 844  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetMatProperties,Err)
 845  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldMatProperties,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 846    & MaterialsFieldUserNumberMatPropertiesPorosity,POROSITY_PARAM_MAT_PROPERTIES,Err)
 847  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldMatProperties,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 848    & MaterialsFieldUserNumberMatPropertiesPermOverVis,PERM_OVER_VIS_PARAM_MAT_PROPERTIES,Err)
 849
 850
 851  !
 852  !================================================================================================================================
 853  !
 854
 855  !INDEPENDENT FIELDS
 856
 857  !Create the equations set independent field variables for the solid
 858  CALL CMISSField_Initialise(IndependentFieldSolid,Err)
 859  CALL CMISSEquationsSet_IndependentCreateStart(EquationsSetSolid,IndependentFieldUserNumberSolid, &
 860    & IndependentFieldSolid,Err)
 861  !Set the mesh component to be used by the field components.
 862  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 863    CALL CMISSField_ComponentMeshComponentSet(IndependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
 864      & DisplacementMeshComponentNumber,Err)
 865  ENDDO
 866  !Finish the equations set independent field variables
 867  CALL CMISSEquationsSet_IndependentCreateFinish(EquationsSetSolid,Err)
 868
 869  !
 870  !================================================================================================================================
 871  !
 872
 873  !EQUATIONS
 874
 875  !Create the equations set equations
 876  CALL CMISSEquations_Initialise(EquationsDarcy,Err)
 877  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy,EquationsDarcy,Err)
 878  !Set the equations matrices sparsity type
 879  CALL CMISSEquations_SparsityTypeSet(EquationsDarcy,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
 880!   !Set the equations lumping type
 881!   CALL CMISSEquations_LumpingTypeSet(EquationsDarcy,CMISS_EQUATIONS_UNLUMPED_MATRICES,Err)
 882  !Set the equations set output
 883  CALL CMISSEquations_OutputTypeSet(EquationsDarcy,EQUATIONS_DARCY_OUTPUT,Err)
 884  !Finish the equations set equations
 885  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetDarcy,Err)
 886
 887  !Create the equations set equations
 888  CALL CMISSEquations_Initialise(EquationsMatProperties,Err)
 889  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetMatProperties,EquationsMatProperties,Err)
 890  !Set the equations matrices sparsity type
 891  CALL CMISSEquations_SparsityTypeSet(EquationsMatProperties,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
 892  !Set the equations set output
 893  CALL CMISSEquations_OutputTypeSet(EquationsMatProperties,EQUATIONS_MAT_PROPERTIES_OUTPUT,Err)
 894  !Finish the equations set equations
 895  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetMatProperties,Err)
 896
 897  !--------------------------------------------------------------------------------------------------------------------------------
 898  ! Solid
 899
 900  !Create the equations set equations
 901  CALL CMISSEquations_Initialise(EquationsSolid,Err)
 902  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetSolid,EquationsSolid,Err)
 903  CALL CMISSEquations_SparsityTypeSet(EquationsSolid,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
 904  CALL CMISSEquations_OutputTypeSet(EquationsSolid,CMISS_EQUATIONS_NO_OUTPUT,Err)
 905  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetSolid,Err)
 906
 907  ! end Solid
 908  !--------------------------------------------------------------------------------------------------------------------------------
 909
 910  !
 911  !================================================================================================================================
 912  !
 913  CALL CMISSField_Initialise(AnalyticFieldDarcy,Err)
 914    CALL CMISSEquationsSet_AnalyticCreateStart(EquationsSetDarcy,&
 915      & CMISS_EQUATIONS_SET_INCOMP_ELAST_DARCY_ANALYTIC_DARCY,&
 916      & AnalyticFieldUserNumberDarcy,AnalyticFieldDarcy,Err)
 917  
 918    !Finish the equations set analytic field variables
 919    CALL CMISSEquationsSet_AnalyticCreateFinish(EquationsSetDarcy,Err)
 920
 921
 922
 923  !--------------------------------------------------------------------------------------------------------------------------------
 924  ! Solid
 925
 926  !Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure
 927  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 928    & 1,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,Err)
 929  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 930    & 2,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,Err)
 931  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
 932    & 3,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,Err)
 933  CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4, &
 934    & 0.0_CMISSDP, &
 935    & Err)
 936
 937  ! end Solid
 938  !--------------------------------------------------------------------------------------------------------------------------------
 939
 940  !
 941  !================================================================================================================================
 942  !
 943
 944  !PROBLEMS
 945
 946  !Start the creation of a problem.
 947  CALL CMISSProblem_Initialise(Problem,Err)
 948  CALL CMISSControlLoop_Initialise(ControlLoop,Err)
 949  CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
 950  !Set the problem to be a ALE Darcy problem
 951  CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_MULTI_PHYSICS_CLASS,CMISS_PROBLEM_FINITE_ELASTICITY_DARCY_TYPE, &
 952    & CMISS_PROBLEM_QUASISTATIC_ELASTICITY_TRANSIENT_DARCY_SUBTYPE,Err)
 953  !Finish the creation of a problem.
 954  CALL CMISSProblem_CreateFinish(Problem,Err)
 955  !Start the creation of the problem control loop
 956  CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
 957  !Get the control loop
 958  CALL CMISSProblem_ControlLoopGet(Problem,CMISS_CONTROL_LOOP_NODE,ControlLoop,Err)
 959!   CALL CMISSControlLoop_MaximumIterationsSet(ControlLoop,2,Err)
 960  !Set the times
 961  CALL CMISSControlLoop_TimesSet(ControlLoop,DYNAMIC_SOLVER_DARCY_START_TIME,DYNAMIC_SOLVER_DARCY_STOP_TIME, &
 962    & DYNAMIC_SOLVER_DARCY_TIME_INCREMENT,Err)
 963  !Set the output timing
 964  CALL CMISSControlLoop_TimeOutputSet(ControlLoop,DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY,Err)
 965!   !Set the output type
 966!   CALL CMISSControlLoop_OutputTypeSet(ControlLoop,CMISS_CONTROL_LOOP_PROGRESS_OUTPUT,Err)
 967  !Finish creating the problem control loop
 968  CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
 969
 970
 971  !
 972  !================================================================================================================================
 973  !
 974
 975  !SOLVERS
 976
 977  !Start the creation of the problem solvers
 978  CALL CMISSSolver_Initialise(SolverSolid,Err)
 979  CALL CMISSSolver_Initialise(LinearSolverMatProperties,Err)
 980  CALL CMISSSolver_Initialise(DynamicSolverDarcy,Err)
 981  CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
 982
 983  CALL CMISSProblem_SolversCreateStart(Problem,Err)
 984
 985  !--------------------------------------------------------------------------------------------------------------------------------
 986  ! Solid
 987
 988  !Get the finite elasticity solver
 989  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
 990    & SolverSolidIndex,SolverSolid,Err)
 991  CALL CMISSSolver_OutputTypeSet(SolverSolid,CMISS_SOLVER_PROGRESS_OUTPUT,Err)
 992!   CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_FD_CALCULATED,Err)
 993  CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_EQUATIONS_CALCULATED,Err)
 994
 995  CALL CMISSSolver_NewtonAbsoluteToleranceSet(SolverSolid,ABSOLUTE_TOLERANCE,Err)
 996  CALL CMISSSolver_NewtonRelativeToleranceSet(SolverSolid,RELATIVE_TOLERANCE,Err)
 997  CALL CMISSSolver_NewtonMaximumIterationsSet(SolverSolid,MAXIMUM_ITERATIONS,Err)
 998
 999!   CALL CMISSSolverNonLinearTypeSet(SolverSolid,CMISS_SOLVER_NONLINEAR_NEWTON,Err)
1000!   CALL CMISSSolver_LibraryTypeSet(SolverSolid,CMISS_SOLVER_PETSC_LIBRARY,Err)
1001
1002!   CALL CMISSSolver_NewtonLinearSolverGet(SolverSolid,LinearSolverSolid,Err)
1003!   CALL CMISSSolver_LinearTypeSet(LinearSolverSolid,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
1004
1005  ! end Solid
1006  !--------------------------------------------------------------------------------------------------------------------------------
1007
1008  !Get the deformation-dependent material properties solver
1009  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
1010    & SolverMatPropertiesIndex,LinearSolverMatProperties,Err)
1011  !Set the output type
1012  CALL CMISSSolver_OutputTypeSet(LinearSolverMatProperties,LINEAR_SOLVER_MAT_PROPERTIES_OUTPUT_TYPE,Err)
1013  !Set the solver settings
1014  IF(LINEAR_SOLVER_MAT_PROPERTIES_DIRECT_FLAG) THEN
1015    CALL CMISSSolver_LinearTypeSet(LinearSolverMatProperties,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
1016    CALL CMISSSolver_LibraryTypeSet(LinearSolverMatProperties,CMISS_SOLVER_MUMPS_LIBRARY,Err)
1017  ELSE
1018    CALL CMISSSolver_LinearTypeSet(LinearSolverMatProperties,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
1019    CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverMatProperties,MAXIMUM_ITERATIONS,Err)
1020    CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverMatProperties,DIVERGENCE_TOLERANCE,Err)
1021    CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverMatProperties,RELATIVE_TOLERANCE,Err)
1022    CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverMatProperties,ABSOLUTE_TOLERANCE,Err)
1023    CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverMatProperties,RESTART_VALUE,Err)
1024  ENDIF
1025
1026  !Get the Darcy solver
1027  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
1028    & SolverDarcyIndex,DynamicSolverDarcy,Err)
1029  !Set the output type
1030  CALL CMISSSolver_OutputTypeSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE,Err)
1031  !Set theta
1032  CALL CMISSSolver_DynamicThetaSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_THETA,Err)
1033!   CALL CMISSSolverDynamicDynamicSet(DynamicSolverDarcy,.TRUE.,Err)
1034  !Get the dynamic linear solver
1035  CALL CMISSSolver_DynamicLinearSolverGet(DynamicSolverDarcy,LinearSolverDarcy,Err)
1036  !Set the solver settings
1037  IF(LINEAR_SOLVER_DARCY_DIRECT_FLAG) THEN
1038    CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
1039    CALL CMISSSolver_LibraryTypeSet(LinearSolverDarcy,CMISS_SOLVER_MUMPS_LIBRARY,Err)
1040  ELSE
1041    CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
1042    CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverDarcy,MAXIMUM_ITERATIONS,Err)
1043    CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverDarcy,DIVERGENCE_TOLERANCE,Err)
1044    CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverDarcy,RELATIVE_TOLERANCE,Err)
1045    CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverDarcy,ABSOLUTE_TOLERANCE,Err)
1046    CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverDarcy,RESTART_VALUE,Err)
1047  ENDIF
1048
1049  !Finish the creation of the problem solver
1050  CALL CMISSProblem_SolversCreateFinish(Problem,Err)
1051
1052  !
1053  !================================================================================================================================
1054  !
1055
1056  !SOLVER EQUATIONS
1057
1058  !Start the creation of the problem solver equations
1059  CALL CMISSSolver_Initialise(SolverSolid,Err)
1060  CALL CMISSSolver_Initialise(LinearSolverMatProperties,Err)
1061  CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
1062
1063  CALL CMISSSolverEquations_Initialise(SolverEquationsSolid,Err)
1064  CALL CMISSSolverEquations_Initialise(SolverEquationsMatProperties,Err)
1065  CALL CMISSSolverEquations_Initialise(SolverEquationsDarcy,Err)
1066
1067  CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
1068  !
1069  !Get the finite elasticity solver equations
1070  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
1071    & SolverSolidIndex,SolverSolid,Err)
1072  CALL CMISSSolver_SolverEquationsGet(SolverSolid,SolverEquationsSolid,Err)
1073  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsSolid,CMISS_SOLVER_SPARSE_MATRICES,Err)
1074  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsSolid,EquationsSetSolid,EquationsSetIndex,Err)
1075  !
1076  !Get the deformation-dependent material properties solver equations
1077  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
1078    & SolverMatPropertiesIndex,LinearSolverMatProperties,Err)
1079  CALL CMISSSolver_SolverEquationsGet(LinearSolverMatProperties,SolverEquationsMatProperties,Err)
1080  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsMatProperties,CMISS_SOLVER_SPARSE_MATRICES,Err)
1081  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsMatProperties,EquationsSetMatProperties,EquationsSetIndex,Err)
1082  !
1083  !Get the Darcy solver equations
1084  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
1085    & SolverDarcyIndex,LinearSolverDarcy,Err)
1086  CALL CMISSSolver_SolverEquationsGet(LinearSolverDarcy,SolverEquationsDarcy,Err)
1087  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsDarcy,CMISS_SOLVER_SPARSE_MATRICES,Err)
1088  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsDarcy,EquationsSetDarcy,EquationsSetIndex,Err)
1089  !
1090  !Finish the creation of the problem solver equations
1091  CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
1092
1093  !
1094  !================================================================================================================================
1095  !
1096
1097  !BOUNDARY CONDITIONS
1098  !Start the creation of the equations set boundary conditions for Darcy
1099  CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsDarcy,Err)
1100  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsDarcy,BoundaryConditionsDarcy,Err)
1101
1102  CALL CMISSSolverEquations_BoundaryConditionsAnalytic(SolverEquationsDarcy,Err)
1103
1104
1105  !Prescribe boundary conditions (absolute nodal parameters)
1106  !Solid is computed in absolute position, rather than displacement. Thus BCs for absolute position
1107  CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsSolid,Err)
1108  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsSolid,BoundaryConditionsSolid,Err)
1109
1110!   !--- BCs on normal velocity only
1111!   CONDITION = CMISS_BOUNDARY_CONDITION_MOVED_WALL
1112! 
1113!   IF( CM%D==2_CMISSIntg ) THEN !CM%D = number of dimensions, ie 2D
1114!     DO NODE_NUMBER=1_CMISSIntg,NUMBER_OF_NODES_GEOMETRY
1115!       COORD_X = CM%N(NODE_NUMBER,1_CMISSIntg)
1116!       COORD_Y = CM%N(NODE_NUMBER,2_CMISSIntg)
1117! 
1118!       IF( (ABS(COORD_X-DOMAIN_X1) < GEOMETRY_TOLERANCE) ) THEN
1119!         !x-velocity
1120!         VALUE = 1.0_CMISSDP
1121!         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1122!           & NODE_NUMBER,1_CMISSIntg,CONDITION,VALUE,Err)
1123!       END IF
1124!       !
1125!       IF( (ABS(COORD_X-DOMAIN_X2) < GEOMETRY_TOLERANCE) ) THEN
1126!         !x-velocity
1127!         VALUE = 1.0_CMISSDP
1128!         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1129!           & NODE_NUMBER,1_CMISSIntg,CONDITION,VALUE,Err)
1130!       END IF
1131!       !
1132!       IF( (ABS(COORD_Y-DOMAIN_Y1) < GEOMETRY_TOLERANCE) ) THEN
1133!         !y-velocity
1134!         VALUE = 2.0_CMISSDP
1135!         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1136!           & NODE_NUMBER,2_CMISSIntg,CONDITION,VALUE,Err)
1137!       END IF
1138!       !
1139!       IF( (ABS(COORD_Y-DOMAIN_Y2) < GEOMETRY_TOLERANCE) ) THEN
1140!         !y-velocity
1141!         VALUE = 2.0_CMISSDP
1142!         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1143!           & NODE_NUMBER,2_CMISSIntg,CONDITION,VALUE,Err)
1144!       END IF
1145!     END DO
1146!   ELSE IF( CM%D==3_CMISSIntg ) THEN ! 3D geometry
1147!     DO NODE_NUMBER=1_CMISSIntg,NUMBER_OF_NODES_GEOMETRY  !What if different number of nodes geometry and velocity ?
1148!       COORD_X = CM%N(NODE_NUMBER,1_CMISSIntg)
1149!       COORD_Y = CM%N(NODE_NUMBER,2_CMISSIntg)
1150!       COORD_Z = CM%N(NODE_NUMBER,3_CMISSIntg)
1151! 
1152!       IF( (ABS(COORD_X-DOMAIN_X1) < GEOMETRY_TOLERANCE) ) THEN
1153! !         !x-velocity: F L U I D ( V Variable type )
1154! !         VALUE = 10.0_CMISSDP
1155! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,,DependentFieldCMISSFieldVVariableType,CMISS_NO_GLOBAL_DERIV, &
1156! ! !           & NODE_NUMBER,4_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err) !BC on pressure component
1157! ! !           & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_MOVED_WALL,VALUE,Err) !inflow
1158! !           & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err) !time-dependent inflow
1159! 
1160! !         !x-position: S O L I D ( U Variable Type)
1161! ! !         VALUE = 1.0_CMISSDP * DOMAIN_X1
1162! !         VALUE = COORD_X
1163! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1164! !           & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1165! 
1166! !         EDGE_COUNT = 0_CMISSIntg
1167! !         IF( (ABS(COORD_Y-DOMAIN_Y1) < GEOMETRY_TOLERANCE) ) EDGE_COUNT = EDGE_COUNT + 1_CMISSIntg
1168! !         IF( (ABS(COORD_Y-DOMAIN_Y2) < GEOMETRY_TOLERANCE) ) EDGE_COUNT = EDGE_COUNT + 1_CMISSIntg
1169! !         IF( (ABS(COORD_Z-DOMAIN_Z1) < GEOMETRY_TOLERANCE) ) EDGE_COUNT = EDGE_COUNT + 1_CMISSIntg
1170! !         IF( (ABS(COORD_Z-DOMAIN_Z2) < GEOMETRY_TOLERANCE) ) EDGE_COUNT = EDGE_COUNT + 1_CMISSIntg
1171! ! 
1172! !         IF(EDGE_COUNT == 2_CMISSIntg) THEN !it is a corner node
1173! !           VALUE = COORD_Y
1174! !           CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1175! !             & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1176! ! 
1177! !           VALUE = COORD_Z
1178! !           CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1179! !             & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1180! !         END IF
1181!       END IF
1182!       !
1183!       IF( (ABS(COORD_X-DOMAIN_X2) < GEOMETRY_TOLERANCE) ) THEN
1184! !         !x-velocity: F L U I D
1185! !         VALUE = 10.0_CMISSDP
1186! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1187! ! !           & NODE_NUMBER,4_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err) !BC on pressure component
1188! ! !           & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_MOVED_WALL,VALUE,Err) !impermeable wall, zero flux
1189! !           & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err) !impermeable wall, zero flux
1190! 
1191! !         !x-position: S O L I D
1192! !         VALUE = COORD_X
1193! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1, &
1194! !           & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1195! !         
1196! !         !Fix point 1
1197! !         IF( (ABS(COORD_Y-DOMAIN_Y2) < GEOMETRY_TOLERANCE) ) THEN
1198! !           IF( (ABS(COORD_Z-DOMAIN_Z2) < GEOMETRY_TOLERANCE) ) THEN
1199! !             VALUE = COORD_Y
1200! !             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1201! !               & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1202! ! 
1203! !             VALUE = COORD_Z
1204! !             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1205! !               & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1206! !           END IF
1207! !         END IF
1208! !         !(Fix) point 2
1209! !         IF( (ABS(COORD_Y-DOMAIN_Y1) < GEOMETRY_TOLERANCE) ) THEN
1210! !           IF( (ABS(COORD_Z-DOMAIN_Z2) < GEOMETRY_TOLERANCE) ) THEN
1211! !             VALUE = COORD_Z
1212! !             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1213! !               & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1214! !           END IF
1215! !         END IF
1216! !         !(Fix) point 3
1217! !         IF( (ABS(COORD_Y-DOMAIN_Y2) < GEOMETRY_TOLERANCE) ) THEN
1218! !           IF( (ABS(COORD_Z-DOMAIN_Z1) < GEOMETRY_TOLERANCE) ) THEN
1219! !             VALUE = COORD_Y
1220! !             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1221! !               & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1222! !           END IF
1223! !         END IF
1224! 
1225! 
1226!       END IF
1227!       !
1228!       IF( (ABS(COORD_Y-DOMAIN_Y1) < GEOMETRY_TOLERANCE) ) THEN
1229! !         !y-velocity: F L U I D
1230! !         VALUE = 0.0_CMISSDP
1231! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1232! !           & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_MOVED_WALL,VALUE,Err)
1233! ! 
1234! !         !y-position: S O L I D
1235! !         VALUE = 1.0_CMISSDP * DOMAIN_Y1
1236! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1237! ! !           & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_MOVED_WALL_INCREMENTED,VALUE,Err)
1238! !           & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1239!       END IF
1240!       !
1241!       IF( (ABS(COORD_Y-DOMAIN_Y2) < GEOMETRY_TOLERANCE) ) THEN
1242! !         !y-velocity: F L U I D
1243! !         VALUE = 0.0_CMISSDP
1244! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1245! !           & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_MOVED_WALL,VALUE,Err)
1246! ! 
1247! ! !         !y-position: S O L I D
1248! ! !         VALUE = 1.1_CMISSDP * DOMAIN_Y2
1249! ! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1, &
1250! ! !           & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1251!       END IF
1252!       !
1253!       IF( (ABS(COORD_Z-DOMAIN_Z1) < GEOMETRY_TOLERANCE) ) THEN
1254!         !z-velocity: F L U I D
1255!         !mass-correction: F L U I D
1256!         VALUE = 0.1_CMISSDP
1257!         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1258! !           & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1259! !           & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_CORRECTION_MASS_INCREASE,VALUE,Err)
1260! !           & NODE_NUMBER,4_CMISSIntg,CMISS_BOUNDARY_CONDITION_CORRECTION_MASS_INCREASE,VALUE,Err)
1261!           & NODE_NUMBER,4_CMISSIntg,CMISS_BOUNDARY_CONDITION_FREE,VALUE,Err)
1262! 
1263! !         !z-position: S O L I D
1264! !         VALUE = 1.0_CMISSDP * DOMAIN_Z1
1265! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1266! !           & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1267!       END IF
1268!       !
1269!       IF( (ABS(COORD_Z-DOMAIN_Z2) < GEOMETRY_TOLERANCE) ) THEN
1270! !         !z-velocity: F L U I D
1271! !         VALUE = 0.0_CMISSDP
1272! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsDarcy,DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1273! !           & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1274! 
1275! !         VALUE = COORD_X
1276! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1, &
1277! !           & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1278! ! 
1279! !         VALUE = COORD_Y
1280! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1281! !           & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1282! ! 
1283! !         VALUE = COORD_Z
1284! !         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1285! !           & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1286! 
1287!         !x-position: S O L I D
1288!         VALUE = COORD_Z
1289!         CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1, &
1290!           & NODE_NUMBER,3_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1291!         
1292!         !Fix point 1
1293!         IF( (ABS(COORD_Y-DOMAIN_Y2) < GEOMETRY_TOLERANCE) ) THEN
1294!           IF( (ABS(COORD_X-DOMAIN_X2) < GEOMETRY_TOLERANCE) ) THEN
1295!             VALUE = COORD_Y
1296!             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1297!               & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1298! 
1299!             VALUE = COORD_X
1300!             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1301!               & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1302!           END IF
1303!         END IF
1304!         !(Fix) point 2
1305!         IF( (ABS(COORD_Y-DOMAIN_Y1) < GEOMETRY_TOLERANCE) ) THEN
1306!           IF( (ABS(COORD_X-DOMAIN_X2) < GEOMETRY_TOLERANCE) ) THEN
1307!             VALUE = COORD_X
1308!             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1309!               & NODE_NUMBER,1_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1310!           END IF
1311!         END IF
1312!         !(Fix) point 3
1313!         IF( (ABS(COORD_Y-DOMAIN_Y2) < GEOMETRY_TOLERANCE) ) THEN
1314!           IF( (ABS(COORD_X-DOMAIN_X1) < GEOMETRY_TOLERANCE) ) THEN
1315!             VALUE = COORD_Y
1316!             CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_NO_GLOBAL_DERIV, &
1317!               & NODE_NUMBER,2_CMISSIntg,CMISS_BOUNDARY_CONDITION_FIXED,VALUE,Err)
1318!           END IF
1319!         END IF
1320! 
1321! 
1322!       END IF
1323!     END DO
1324!   END IF
1325
1326  !Finish the creation of the equations set boundary conditions for Darcy
1327  !CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsDarcy,Err)
1328  !Finish the creation of the equations set boundary conditions for the solid
1329  CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsSolid,Err)
1330  !
1331  !Start the creation of the equations set boundary conditions for deformation-dependent material properties
1332  CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsMatProperties,Err)
1333  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsMatProperties,BoundaryConditionsMatProperties,Err)
1334  !(No boundary conditions requrired for deformation-dependent material properties)
1335  !Finish the creation of the equations set boundary conditions for deformation-dependent material properties
1336  CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsMatProperties,Err)
1337
1338
1339  !
1340  !================================================================================================================================
1341  !
1342
1343  !RUN SOLVERS
1344
1345  !Turn of PETSc error handling
1346  !CALL PETSC_ERRORHANDLING_SET_ON(ERR,ERROR,*999)
1347
1348  !Solve the problem
1349  WRITE(*,'(A)') "Solving problem..."
1350  CALL CMISSProblem_Solve(Problem,Err)
1351  WRITE(*,'(A)') "Problem solved!"
1352
1353
1354  !
1355  !================================================================================================================================
1356  !
1357
1358  !OUTPUT
1359
1360  EXPORT_FIELD_IO=.FALSE.
1361  IF(EXPORT_FIELD_IO) THEN
1362    WRITE(*,'(A)') "Exporting fields..."
1363    CALL CMISSFields_Initialise(Fields,Err)
1364    CALL CMISSFields_Create(Region,Fields,Err)
1365    CALL CMISSFields_NodesExport(Fields,"FiniteElasticityDarcy","FORTRAN",Err)
1366    CALL CMISSFields_ElementsExport(Fields,"FiniteElasticityDarcy","FORTRAN",Err)
1367    CALL CMISSFields_Finalise(Fields,Err)
1368    WRITE(*,'(A)') "Field exported!"
1369  ENDIF
1370
1371
1372  !Finialise CMISS
1373!   CALL CMISSFinalise(Err)
1374
1375  WRITE(*,'(A)') "Program successfully completed."
1376
1377  STOP
1378
1379END PROGRAM INCOMPELASTDRIVENDARCYANALYTICDARCYEXAMPLE