PageRenderTime 92ms CodeModel.GetById 23ms app.highlight 60ms RepoModel.GetById 1ms app.codeStats 1ms

/MultiPhysics/CoupledDiffusionDiffusion/MonolithicSchemeTestAnalytic/src/MonolithicSchemeTestAnalyticExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 1129 lines | 516 code | 127 blank | 486 comment | 0 complexity | b0aa9e05b74e8f26ad63c016923ddb23 MD5 | raw file
   1!> \file
   2!> \authors Andrew Cookson
   3!> \brief This is an example program to solve coupled multi-compartment diffusion equations in monolithic scheme 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/CoupledDiffusionDiffusion/MonolithicSchemeTest/src/MonolithicSchemeTestExample.f90
  43!! Example program to solve coupled MonolithicSchemeTest 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/CoupledDiffusionDiffusion/MonolithicSchemeTest/build-intel'>Linux Intel Build</a>
  46!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/MultiPhysics/CoupledDiffusionDiffusion/MonolithicSchemeTest/build-intel'>Linux GNU Build</a>
  47!!
  48!<
  49
  50! ! 
  51! !  This example considers a volume coupled multi-compartment diffusion model, as a means of testing the monolithic assembly of
  52! !  such a system of equations, for use in more complicated problem in future.
  53! !  This example will initially couple together three diffusion equations, with the transfer between each equation being proportional to
  54! !  the concentration difference between each equation.
  55
  56!> Main program
  57
  58PROGRAM MONOLITHICSCHEMETESTFIELDMLEXAMPLE
  59
  60  !
  61  !================================================================================================================================
  62  !
  63
  64  !PROGRAM LIBRARIES
  65
  66  USE OPENCMISS
  67  USE FIELDML_OUTPUT_ROUTINES
  68  USE FIELDML_UTIL_ROUTINES
  69  USE FIELDML_API
  70  USE FLUID_MECHANICS_IO_ROUTINES
  71  USE MPI
  72
  73#ifdef WIN32
  74  USE IFQWINCMISS
  75#endif
  76
  77  !
  78  !================================================================================================================================
  79  !
  80
  81  !PROGRAM VARIABLES AND TYPES
  82
  83  IMPLICIT NONE
  84
  85  !Test program parameters
  86
  87  REAL(CMISSDP), PARAMETER :: HEIGHT=1.0_CMISSDP
  88  REAL(CMISSDP), PARAMETER :: WIDTH=1.0_CMISSDP
  89  REAL(CMISSDP), PARAMETER :: LENGTH=3.0_CMISSDP
  90
  91  INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
  92  INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=2
  93  INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=3
  94  INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=4
  95  INTEGER(CMISSIntg), PARAMETER :: GeometricFieldUserNumber=5
  96  INTEGER(CMISSIntg), PARAMETER :: DependentFieldUserNumber=6
  97  INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=9
  98  INTEGER(CMISSIntg) :: MaterialsFieldUserNumberDiffusion
  99!   INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDiffusionTwo=10
 100!   INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDiffusionThree=11
 101  INTEGER(CMISSIntg) :: SourceFieldUserNumberDiffusion
 102!   INTEGER(CMISSIntg), PARAMETER :: SourceFieldUserNumberDiffusionTwo=16
 103!   INTEGER(CMISSIntg), PARAMETER :: SourceFieldUserNumberDiffusionThree=17
 104  INTEGER(CMISSIntg) :: EquationsSetUserNumberDiffusion
 105  INTEGER(CMISSIntg) :: AnalyticFieldUserNumberDiffusion
 106!   INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberDiffusionTwo=13
 107!   INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberDiffusionThree=14
 108  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=20
 109!!!!!!!!!!!!!!!!!!!EQUATIONS FIELD PARAMETERS
 110  INTEGER(CMISSIntg), PARAMETER :: EquationsFieldDiffusionOne=33
 111  INTEGER(CMISSIntg), PARAMETER :: EquationsFieldDiffusionTwo=34
 112  INTEGER(CMISSIntg), PARAMETER :: EquationsFieldDiffusionThree=35
 113
 114
 115
 116
 117
 118  INTEGER(CMISSIntg) :: EquationsSetFieldUserNumberDiffusion
 119  INTEGER(CMISSIntg) :: icompartment,Ncompartments,num_var
 120!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 121  INTEGER(CMISSIntg), PARAMETER :: DomainUserNumber=1
 122  INTEGER(CMISSIntg), PARAMETER :: NumberOfUserDomains=1
 123  INTEGER(CMISSIntg), PARAMETER :: SolverDiffusionUserNumber=1
 124  !Program types
 125
 126  TYPE(EXPORT_CONTAINER):: CM
 127
 128  !Program variables
 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) :: NUMBER_OF_COMPARTMENTS
 137  
 138  INTEGER(CMISSIntg) :: BASIS_TYPE
 139  INTEGER(CMISSIntg) :: BASIS_NUMBER_GEOMETRY
 140  INTEGER(CMISSIntg) :: BASIS_NUMBER_CONC_ONE
 141  INTEGER(CMISSIntg) :: BASIS_NUMBER_CONC_TWO
 142  INTEGER(CMISSIntg) :: BASIS_NUMBER_CONC_THREE
 143  INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_GEOMETRY
 144  INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_CONC_ONE
 145  INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_CONC_TWO
 146  INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_CONC_THREE
 147  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_GEOMETRY
 148  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_CONC_ONE
 149  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_CONC_TWO
 150  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_CONC_THREE
 151  INTEGER(CMISSIntg) :: MESH_NUMBER_OF_COMPONENTS,MESH_NUMBER_OF_ALL_COMPONENTS
 152  INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_GEOMETRY
 153  INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_CONC_ONE
 154  INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_CONC_TWO
 155  INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_CONC_THREE
 156  INTEGER(CMISSIntg) :: NUMBER_OF_NODES_GEOMETRY
 157  INTEGER(CMISSIntg) :: NUMBER_OF_NODES_CONC_ONE
 158  INTEGER(CMISSIntg) :: NUMBER_OF_NODES_CONC_TWO
 159  INTEGER(CMISSIntg) :: NUMBER_OF_NODES_CONC_THREE
 160  INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_GEOMETRY
 161  INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_CONC_ONE
 162  INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_CONC_TWO
 163  INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_CONC_THREE
 164  INTEGER(CMISSIntg) :: TOTAL_NUMBER_OF_NODES,TOTAL_NUMBER_OF_ALL_NODES
 165  INTEGER(CMISSIntg) :: TOTAL_NUMBER_OF_ELEMENTS
 166  INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
 167  INTEGER(CMISSIntg) :: RESTART_VALUE
 168  INTEGER(CMISSIntg) :: NUMBER_OF_FIXED_WALL_NODES_DIFFUSION_ONE
 169  INTEGER(CMISSIntg) :: NUMBER_OF_INLET_WALL_NODES_DIFFUSION_ONE
 170  INTEGER(CMISSIntg) :: NUMBER_OF_FIXED_WALL_NODES_DIFFUSION_TWO
 171  INTEGER(CMISSIntg) :: NUMBER_OF_INLET_WALL_NODES_DIFFUSION_TWO
 172  INTEGER(CMISSIntg) :: NUMBER_OF_FIXED_WALL_NODES_DIFFUSION_THREE
 173  INTEGER(CMISSIntg) :: NUMBER_OF_INLET_WALL_NODES_DIFFUSION_THREE
 174  INTEGER(CMISSIntg) :: EQUATIONS_DIFFUSION_OUTPUT
 175  INTEGER(CMISSIntg) :: EQUATIONS_DIFFUSION_TWO_OUTPUT
 176  INTEGER(CMISSIntg) :: EQUATIONS_DIFFUSION_THREE_OUTPUT
 177  INTEGER(CMISSIntg) :: COMPONENT_NUMBER
 178  INTEGER(CMISSIntg) :: NODE_NUMBER
 179  INTEGER(CMISSIntg) :: ELEMENT_NUMBER
 180  INTEGER(CMISSIntg) :: NODE_COUNTER
 181  INTEGER(CMISSIntg) :: CONDITION
 182
 183  INTEGER(CMISSIntg) :: LINEAR_SOLVER_DIFFUSION_OUTPUT_FREQUENCY
 184  INTEGER(CMISSIntg) :: LINEAR_SOLVER_DIFFUSION_OUTPUT_TYPE
 185
 186
 187  REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2
 188  REAL(CMISSDP) :: GEOMETRY_TOLERANCE
 189
 190  INTEGER, ALLOCATABLE, DIMENSION(:):: FIXED_WALL_NODES_DIFFUSION_ONE
 191  INTEGER, ALLOCATABLE, DIMENSION(:):: INLET_WALL_NODES_DIFFUSION_ONE
 192  INTEGER, ALLOCATABLE, DIMENSION(:):: FIXED_WALL_NODES_DIFFUSION_TWO
 193  INTEGER, ALLOCATABLE, DIMENSION(:):: INLET_WALL_NODES_DIFFUSION_TWO
 194  INTEGER, ALLOCATABLE, DIMENSION(:):: FIXED_WALL_NODES_DIFFUSION_THREE
 195  INTEGER, ALLOCATABLE, DIMENSION(:):: INLET_WALL_NODES_DIFFUSION_THREE
 196
 197  REAL(CMISSDP) :: INITIAL_FIELD_DIFFUSION_ONE
 198  REAL(CMISSDP) :: INITIAL_FIELD_DIFFUSION_TWO
 199  REAL(CMISSDP) :: INITIAL_FIELD_DIFFUSION_THREE
 200  REAL(CMISSDP) :: BOUNDARY_CONDITIONS_DIFFUSION_ONE
 201  REAL(CMISSDP) :: BOUNDARY_CONDITIONS_DIFFUSION_TWO
 202  REAL(CMISSDP) :: BOUNDARY_CONDITIONS_DIFFUSION_THREE
 203  REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
 204  REAL(CMISSDP) :: RELATIVE_TOLERANCE
 205  REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
 206  REAL(CMISSDP) :: LINESEARCH_ALPHA
 207  REAL(CMISSDP) :: VALUE
 208
 209  REAL(CMISSDP) :: LINEAR_SOLVER_DIFFUSION_START_TIME
 210  REAL(CMISSDP) :: LINEAR_SOLVER_DIFFUSION_STOP_TIME
 211  REAL(CMISSDP) :: LINEAR_SOLVER_DIFFUSION_TIME_INCREMENT
 212
 213  LOGICAL :: EXPORT_FIELD_IO
 214  LOGICAL :: LINEAR_SOLVER_DIFFUSION_DIRECT_FLAG
 215  LOGICAL :: INLET_WALL_NODES_DIFFUSION_ONE_FLAG
 216  LOGICAL :: INLET_WALL_NODES_DIFFUSION_TWO_FLAG
 217  LOGICAL :: INLET_WALL_NODES_DIFFUSION_THREE_FLAG
 218  !CMISS variables
 219  CHARACTER(C_CHAR), PARAMETER :: NUL=C_NULL_CHAR
 220  !Regions
 221  TYPE(CMISSRegionType) :: Region
 222  TYPE(CMISSRegionType) :: WorldRegion
 223  !Coordinate systems
 224  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem
 225  TYPE(CMISSCoordinateSystemType) :: WorldCoordinateSystem
 226  !Basis
 227  TYPE(CMISSBasisType) :: BasisGeometry
 228  TYPE(CMISSBasisType) :: BasisConcOne
 229  TYPE(CMISSBasisType) :: BasisConcTwo
 230  TYPE(CMISSBasisType) :: BasisConcThree
 231  
 232  !TYPE(CMISSBasisType), ALLOCATABLE, DIMENSION(:) :: BasicConc
 233
 234  !Nodes
 235  TYPE(CMISSNodesType) :: Nodes
 236  !Elements
 237  TYPE(CMISSMeshElementsType) :: MeshElementsGeometry
 238  TYPE(CMISSMeshElementsType) :: MeshElementsConcOne
 239  TYPE(CMISSMeshElementsType) :: MeshElementsConcTwo
 240  TYPE(CMISSMeshElementsType) :: MeshElementsConcThree
 241
 242  TYPE(CMISSMeshElementsType), ALLOCATABLE, DIMENSION(:) :: MeshElementsConc
 243  !Meshes
 244  TYPE(CMISSMeshType) :: Mesh
 245  !Decompositions
 246  TYPE(CMISSDecompositionType) :: Decomposition
 247  !Fields
 248  TYPE(CMISSFieldsType) :: Fields
 249  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh  
 250  !Field types
 251  TYPE(CMISSFieldType) :: GeometricField
 252  TYPE(CMISSFieldType) :: DependentField
 253  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: MaterialsFieldDiffusion
 254  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: SourceFieldDiffusion
 255  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: AnalyticFieldDiffusion
 256  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDiffusion
 257  TYPE(CMISSEquationsSetType), ALLOCATABLE, DIMENSION(:) :: EquationsSetDiffusion
 258  TYPE(CMISSEquationsType), ALLOCATABLE, DIMENSION(:) :: EquationsDiffusion
 259  TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: EquationsSetFieldDiffusion
 260  !Problems
 261  TYPE(CMISSProblemType) :: Problem
 262  !Control loops
 263  TYPE(CMISSControlLoopType) :: ControlLoop
 264  !Solvers
 265  TYPE(CMISSSolverType) :: SolverDiffusion
 266  TYPE(CMISSSolverType) :: LinearSolverDiffusion
 267  !Solver equations
 268  TYPE(CMISSSolverEquationsType) :: SolverEquationsDiffusion
 269
 270#ifdef WIN32
 271  !Quickwin type
 272  LOGICAL :: QUICKWIN_STATUS=.FALSE.
 273  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
 274#endif
 275  
 276  !Generic CMISS variables
 277  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber
 278  INTEGER(CMISSIntg) :: EquationsSetIndex
 279  INTEGER(CMISSIntg) :: EquationsSetIndexTwo
 280  INTEGER(CMISSIntg) :: EquationsSetIndexThree
 281  INTEGER(CMISSIntg) :: Err
 282
 283
 284  !Array containing the field variable types that will be used (for ease of incorporating inside a loop)
 285  INTEGER(CMISSIntg), ALLOCATABLE, DIMENSION(:) :: VariableTypes
 286  REAL(CMISSDP), ALLOCATABLE, DIMENSION(:,:) :: CouplingCoeffs
 287
 288  INTEGER(CMISSIntg) :: DIAG_LEVEL_LIST(5)
 289!   CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(8) !,TIMING_ROUTINE_LIST(1)
 290  CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(1) !,TIMING_ROUTINE_LIST(1)
 291 
 292  INTEGER(CMISSIntg) :: TotalNumberOfSolidNodes
 293!   INTEGER(CMISSIntg) :: NumberOfSolidMeshComponents
 294  !
 295  !--------------------------------------------------------------------------------------------------------------------------------
 296  !
 297
 298
 299  !FieldML variables
 300  CHARACTER(KIND=C_CHAR,LEN=*), PARAMETER :: outputDirectory = ""
 301  CHARACTER(KIND=C_CHAR,LEN=*), PARAMETER :: outputFilename = "MonolithicMultiCompDiffusion.xml"
 302  CHARACTER(KIND=C_CHAR,LEN=*), PARAMETER :: basename = "monolithic_multicomp_diffusion"
 303
 304  TYPE(FieldmlInfoType) :: fieldmlInfo
 305
 306
 307
 308
 309#ifdef WIN32
 310  !Initialise QuickWin
 311  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
 312  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
 313  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
 314  !Set the window parameters
 315  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
 316  !If attempt fails set with system estimated values
 317  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
 318#endif
 319
 320  !
 321  !================================================================================================================================
 322  !
 323  !INITIALISE OPENCMISS
 324  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
 325  !CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
 326  !
 327  !================================================================================================================================
 328  !
 329  !PROBLEM CONTROL PANEL
 330  !Import cmHeart mesh information
 331  CALL FLUID_MECHANICS_IO_READ_CMHEART(CM,Err)  
 332  BASIS_NUMBER_GEOMETRY=CM%ID_M
 333  BASIS_NUMBER_CONC_ONE=CM%ID_V !USE THE V COMPONENT OF CMHEART INPUT FOR CONCENTRATION ONE
 334  BASIS_NUMBER_CONC_TWO=CM%ID_P !USE THE p COMPONENT OF CMHEART INPUT FOR CONCENTRATION TWO
 335  NUMBER_OF_DIMENSIONS=CM%D
 336  BASIS_TYPE=CM%IT_T
 337  BASIS_XI_INTERPOLATION_GEOMETRY=CM%IT_M
 338!   BASIS_XI_INTERPOLATION_GEOMETRY=2
 339  BASIS_XI_INTERPOLATION_CONC_ONE=CM%IT_V
 340  BASIS_XI_INTERPOLATION_CONC_TWO=CM%IT_P
 341  NUMBER_OF_NODES_GEOMETRY=CM%N_M
 342  NUMBER_OF_NODES_CONC_ONE=CM%N_V
 343  NUMBER_OF_NODES_CONC_TWO=CM%N_P
 344  TOTAL_NUMBER_OF_NODES=CM%N_T
 345  TOTAL_NUMBER_OF_ELEMENTS=CM%E_T
 346  NUMBER_OF_ELEMENT_NODES_GEOMETRY=CM%EN_M
 347  NUMBER_OF_ELEMENT_NODES_CONC_ONE=CM%EN_V
 348  NUMBER_OF_ELEMENT_NODES_CONC_TWO=CM%EN_P
 349!   !Set domain dimensions
 350!   DOMAIN_X1 =  0.0_CMISSDP
 351!   DOMAIN_X2 =  1.0_CMISSDP
 352!   DOMAIN_Y1 =  0.0_CMISSDP
 353!   DOMAIN_Y2 =  1.0_CMISSDP
 354!   DOMAIN_Z1 =  0.0_CMISSDP
 355!   DOMAIN_Z2 =  1.0_CMISSDP
 356!   !Set geometric tolerance
 357!   GEOMETRY_TOLERANCE = 1.0E-12_CMISSDP
 358!   !Set initial values
 359!   INITIAL_FIELD_DIFFUSION_ONE=1.0_CMISSDP
 360!   INITIAL_FIELD_DIFFUSION_TWO=1.0_CMISSDP
 361!   INITIAL_FIELD_DIFFUSION_THREE=1.0_CMISSDP
 362!   !Set initial boundary conditions
 363!   INLET_WALL_NODES_DIFFUSION_ONE_FLAG=.TRUE.
 364!   IF(INLET_WALL_NODES_DIFFUSION_ONE_FLAG) THEN
 365!     NUMBER_OF_INLET_WALL_NODES_DIFFUSION_ONE=36
 366!     ALLOCATE(INLET_WALL_NODES_DIFFUSION_ONE(NUMBER_OF_INLET_WALL_NODES_DIFFUSION_ONE))
 367!     INLET_WALL_NODES_DIFFUSION_ONE=(/191,155,119,83,23,21,192,156,120,84,24,22,&
 368!      & 198,162,126,90,36,35,204,168,132,96,48,47,210,174,138,102,60,59,216,180,144,108,72,71/)
 369!     !Set initial boundary conditions
 370!     BOUNDARY_CONDITIONS_DIFFUSION_ONE=1.0_CMISSDP
 371!   ENDIF  !Set material parameters
 372! 
 373!   INLET_WALL_NODES_DIFFUSION_TWO_FLAG=.TRUE.
 374!   IF(INLET_WALL_NODES_DIFFUSION_TWO_FLAG) THEN
 375!     NUMBER_OF_INLET_WALL_NODES_DIFFUSION_TWO=36
 376!     ALLOCATE(INLET_WALL_NODES_DIFFUSION_TWO(NUMBER_OF_INLET_WALL_NODES_DIFFUSION_TWO))
 377!     INLET_WALL_NODES_DIFFUSION_TWO=(/191,155,119,83,23,21,192,156,120,84,24,22,&
 378!      & 198,162,126,90,36,35,204,168,132,96,48,47,210,174,138,102,60,59,216,180,144,108,72,71/)
 379!     !Set initial boundary conditions
 380!     BOUNDARY_CONDITIONS_DIFFUSION_TWO=1.0_CMISSDP
 381!   ENDIF  !Set material parameters
 382! 
 383!   INLET_WALL_NODES_DIFFUSION_THREE_FLAG=.TRUE.
 384!   IF(INLET_WALL_NODES_DIFFUSION_THREE_FLAG) THEN
 385!     NUMBER_OF_INLET_WALL_NODES_DIFFUSION_THREE=36
 386!     ALLOCATE(INLET_WALL_NODES_DIFFUSION_THREE(NUMBER_OF_INLET_WALL_NODES_DIFFUSION_THREE))
 387!     INLET_WALL_NODES_DIFFUSION_THREE=(/191,155,119,83,23,21,192,156,120,84,24,22,&
 388!      & 198,162,126,90,36,35,204,168,132,96,48,47,210,174,138,102,60,59,216,180,144,108,72,71/)
 389!     !Set initial boundary conditions
 390!     BOUNDARY_CONDITIONS_DIFFUSION_THREE=1.0_CMISSDP
 391!   ENDIF  !Set material parameters
 392
 393  !Get the computational nodes information
 394  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
 395  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
 396
 397  NUMBER_GLOBAL_X_ELEMENTS=40
 398  NUMBER_GLOBAL_Y_ELEMENTS=40
 399  NUMBER_GLOBAL_Z_ELEMENTS=0
 400
 401  IF(NUMBER_GLOBAL_Z_ELEMENTS==0)THEN
 402    NUMBER_OF_DIMENSIONS=2
 403  ELSE
 404    NUMBER_OF_DIMENSIONS=3
 405  ENDIF
 406  NUMBER_OF_DOMAINS=NumberOfComputationalNodes
 407
 408  !Set all diganostic levels on for testing
 409
 410  CALL MPI_BCAST(NUMBER_GLOBAL_X_ELEMENTS,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
 411  CALL MPI_BCAST(NUMBER_GLOBAL_Y_ELEMENTS,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
 412  CALL MPI_BCAST(NUMBER_GLOBAL_Z_ELEMENTS,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
 413  CALL MPI_BCAST(NUMBER_OF_DOMAINS,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
 414
 415
 416
 417
 418  !Set material parameters
 419  !Set number of Gauss points (Mind that also material field may be interpolated)
 420  BASIS_XI_GAUSS_GEOMETRY=3 !4
 421  BASIS_XI_GAUSS_CONC_ONE=3 !4
 422  BASIS_XI_GAUSS_CONC_TWO=3 !4
 423  !Set output parameter
 424  !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
 425  LINEAR_SOLVER_DIFFUSION_OUTPUT_TYPE=CMISS_SOLVER_MATRIX_OUTPUT
 426  !LINEAR_SOLVER_DIFFUSION_OUTPUT_TYPE=CMISS_SOLVER_NO_OUTPUT
 427  !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
 428  EQUATIONS_DIFFUSION_OUTPUT=CMISS_EQUATIONS_MATRIX_OUTPUT
 429  !EQUATIONS_DIFFUSION_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
 430  EQUATIONS_DIFFUSION_TWO_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
 431  EQUATIONS_DIFFUSION_THREE_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
 432  !Set time parameter
 433  LINEAR_SOLVER_DIFFUSION_START_TIME=0.0_CMISSDP
 434  LINEAR_SOLVER_DIFFUSION_STOP_TIME=0.3000001_CMISSDP 
 435  LINEAR_SOLVER_DIFFUSION_TIME_INCREMENT=0.005_CMISSDP
 436  !Set result output parameter
 437  LINEAR_SOLVER_DIFFUSION_OUTPUT_FREQUENCY=1
 438  !Set solver parameters
 439  LINEAR_SOLVER_DIFFUSION_DIRECT_FLAG=.FALSE.
 440
 441  RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
 442  ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
 443  DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
 444  MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
 445  RESTART_VALUE=30_CMISSIntg !default: 30
 446  LINESEARCH_ALPHA=1.0_CMISSDP
 447
 448
 449  icompartment =1_CMISSIntg
 450  Ncompartments=2_CMISSIntg
 451  !
 452  !================================================================================================================================
 453  !
 454
 455  !Set diagnostics
 456
 457   !DIAG_LEVEL_LIST(1)=1
 458!   DIAG_LEVEL_LIST(2)=2
 459!   DIAG_LEVEL_LIST(3)=3
 460!   DIAG_LEVEL_LIST(4)=4
 461!   DIAG_LEVEL_LIST(5)=5
 462
 463!   DIAG_ROUTINE_LIST(1)="DIFFUSION_EQUATION_FINITE_ELEMENT_CALCULATE"
 464!   DIAG_ROUTINE_LIST(2)="PROBLEM_SOLVER_EQUATIONS_SOLVE"
 465   !DIAG_ROUTINE_LIST(1)="SOLVER_DYNAMIC_CREATE_FINISH"
 466
 467  !CMISS_ALL_DIAG_TYPE/CMISS_IN_DIAG_TYPE/CMISS_FROM_DIAG_TYPE
 468   !CALL CMISSDiagnosticsSetOn(CMISS_IN_DIAG_TYPE,DIAG_LEVEL_LIST,"Diagnostics",DIAG_ROUTINE_LIST,Err)
 469
 470  !CMISS_ALL_TIMING_TYPE/CMISS_IN_TIMING_TYPE/CMISS_FROM_TIMING_TYPE
 471  !TIMING_ROUTINE_LIST(1)="PROBLEM_FINITE_ELEMENT_CALCULATE"
 472  !CALL TIMING_SET_ON(IN_TIMING_TYPE,.TRUE.,"",TIMING_ROUTINE_LIST,ERR,ERROR,*999)
 473
 474  !ALLOCATE THE ARRAYS
 475  ALLOCATE (EquationsSetDiffusion(Ncompartments))
 476  ALLOCATE (EquationsSetFieldDiffusion(Ncompartments))
 477  ALLOCATE (MaterialsFieldDiffusion(Ncompartments))
 478  ALLOCATE (SourceFieldDiffusion(Ncompartments))
 479  ALLOCATE (EquationsDiffusion(Ncompartments))
 480  ALLOCATE (AnalyticFieldDiffusion(Ncompartments))
 481
 482!
 483  !================================================================================================================================
 484  !
 485  !COORDINATE SYSTEM
 486  !Start the creation of a new RC coordinate system
 487  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
 488  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
 489  !Set the coordinate system dimension
 490  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NUMBER_OF_DIMENSIONS,Err)
 491  !Finish the creation of the coordinate system
 492  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
 493  !
 494  !================================================================================================================================
 495  !
 496  !REGION
 497  !For a volume-coupled problem, both concentrations are based in the same region
 498  !Start the creation of a new region
 499  CALL CMISSRegion_Initialise(Region,Err)
 500  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
 501  !Set the regions coordinate system as defined above
 502  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
 503  !Finish the creation of the region
 504  CALL CMISSRegion_CreateFinish(Region,Err)
 505  !
 506  !================================================================================================================================
 507  !
 508  !BASES
 509  !Start the creation of new bases: Geometry
 510  MESH_NUMBER_OF_COMPONENTS=1
 511  CALL CMISSBasis_Initialise(BasisGeometry,Err)
 512  CALL CMISSBasis_CreateStart(BASIS_NUMBER_GEOMETRY,BasisGeometry,Err)
 513  !Set the basis type (Lagrange/Simplex)
 514  CALL CMISSBasis_TypeSet(BasisGeometry,BASIS_TYPE,Err)
 515  !Set the basis xi number
 516  CALL CMISSBasis_NumberOfXiSet(BasisGeometry,NUMBER_OF_DIMENSIONS,Err)
 517  !Set the basis xi interpolation and number of Gauss points
 518  IF(NUMBER_OF_DIMENSIONS==2) THEN
 519    CALL CMISSBasis_InterpolationXiSet(BasisGeometry,(/BASIS_XI_INTERPOLATION_GEOMETRY,BASIS_XI_INTERPOLATION_GEOMETRY/),Err)
 520    CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisGeometry,(/BASIS_XI_GAUSS_GEOMETRY,BASIS_XI_GAUSS_GEOMETRY/),Err)
 521  ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
 522    CALL CMISSBasis_InterpolationXiSet(BasisGeometry,(/BASIS_XI_INTERPOLATION_GEOMETRY,BASIS_XI_INTERPOLATION_GEOMETRY, & 
 523      & BASIS_XI_INTERPOLATION_GEOMETRY/),Err)                         
 524    CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisGeometry,(/BASIS_XI_GAUSS_GEOMETRY,BASIS_XI_GAUSS_GEOMETRY, &
 525      & BASIS_XI_GAUSS_GEOMETRY/),Err)
 526  ENDIF
 527  !Finish the creation of the basis
 528  CALL CMISSBasis_CreateFinish(BasisGeometry,Err)
 529  !
 530  !Start the creation of another basis: Concentration_One
 531!  IF(BASIS_XI_INTERPOLATION_CONC_ONE==BASIS_XI_INTERPOLATION_GEOMETRY) THEN
 532    BasisConcOne=BasisGeometry
 533!  ELSE
 534!     MESH_NUMBER_OF_COMPONENTS=MESH_NUMBER_OF_COMPONENTS+1
 535!     !Initialise a new velocity basis
 536!     CALL CMISSBasis_Initialise(BasisConcOne,Err)
 537!     !Start the creation of a basis
 538!     CALL CMISSBasis_CreateStart(BASIS_NUMBER_CONC_ONE,BasisConcOne,Err)
 539!     !Set the basis type (Lagrange/Simplex)
 540!     CALL CMISSBasis_TypeSet(BasisConcOne,BASIS_TYPE,Err)
 541!     !Set the basis xi number
 542!     CALL CMISSBasis_NumberOfXiSet(BasisConcOne,NUMBER_OF_DIMENSIONS,Err)
 543!     !Set the basis xi interpolation and number of Gauss points
 544!     IF(NUMBER_OF_DIMENSIONS==2) THEN
 545!       CALL CMISSBasis_InterpolationXiSet(BasisConcOne,(/BASIS_XI_INTERPOLATION_CONC_ONE,BASIS_XI_INTERPOLATION_CONC_ONE/),Err)
 546!       CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisConcOne,(/BASIS_XI_GAUSS_CONC_ONE,BASIS_XI_GAUSS_CONC_ONE/),Err)
 547!     ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
 548!       CALL CMISSBasis_InterpolationXiSet(BasisConcOne,(/BASIS_XI_INTERPOLATION_CONC_ONE,BASIS_XI_INTERPOLATION_CONC_ONE, & 
 549!         & BASIS_XI_INTERPOLATION_CONC_ONE/),Err)                         
 550!       CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisConcOne,(/BASIS_XI_GAUSS_CONC_ONE,BASIS_XI_GAUSS_CONC_ONE, & 
 551!         & BASIS_XI_GAUSS_CONC_ONE/),Err)
 552!     ENDIF
 553!     !Finish the creation of the basis
 554!     CALL CMISSBasis_CreateFinish(BasisConcOne,Err)
 555!   ENDIF
 556  !
 557  !Start the creation of another basis: Concentration_Two
 558!  IF(BASIS_XI_INTERPOLATION_CONC_TWO==BASIS_XI_INTERPOLATION_GEOMETRY) THEN
 559    BasisConcTwo=BasisGeometry
 560!  ELSE IF(BASIS_XI_INTERPOLATION_CONC_TWO==BASIS_XI_INTERPOLATION_CONC_ONE) THEN
 561!     BasisConcTwo=BasisConcOne
 562!   ELSE
 563!     MESH_NUMBER_OF_COMPONENTS=MESH_NUMBER_OF_COMPONENTS+1
 564!     !Initialise a new concentration basis
 565!     CALL CMISSBasis_Initialise(BasisConcTwo,Err)
 566!     !Start the creation of a basis
 567!     CALL CMISSBasis_CreateStart(BASIS_NUMBER_CONC_TWO,BasisConcTwo,Err)
 568!     !Set the basis type (Lagrange/Simplex)
 569!     CALL CMISSBasis_TypeSet(BasisConcTwo,BASIS_TYPE,Err)
 570!     !Set the basis xi number
 571!     CALL CMISSBasis_NumberOfXiSet(BasisConcTwo,NUMBER_OF_DIMENSIONS,Err)
 572!     !Set the basis xi interpolation and number of Gauss points
 573!     IF(NUMBER_OF_DIMENSIONS==2) THEN
 574!       CALL CMISSBasis_InterpolationXiSet(BasisConcTwo,(/BASIS_XI_INTERPOLATION_CONC_TWO,BASIS_XI_INTERPOLATION_CONC_TWO/),Err)
 575!       CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisConcTwo,(/BASIS_XI_GAUSS_CONC_TWO,BASIS_XI_GAUSS_CONC_TWO/),Err)
 576!     ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
 577!       CALL CMISSBasis_InterpolationXiSet(BasisConcTwo,(/BASIS_XI_INTERPOLATION_CONC_TWO,BASIS_XI_INTERPOLATION_CONC_TWO, & 
 578!         & BASIS_XI_INTERPOLATION_CONC_TWO/),Err)                         
 579!       CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisConcTwo,(/BASIS_XI_GAUSS_CONC_TWO,BASIS_XI_GAUSS_CONC_TWO, & 
 580!         & BASIS_XI_GAUSS_CONC_TWO/),Err)
 581!     ENDIF
 582!     !Finish the creation of the basis
 583!     CALL CMISSBasis_CreateFinish(BasisConcTwo,Err)
 584!   ENDIF
 585
 586    BasisConcThree=BasisGeometry
 587  !
 588  !================================================================================================================================
 589  !
 590!   !MESH
 591!   !All types of physics utilize the same "mesh", but may be represented on individual mesh components.
 592!   TotalNumberOfSolidNodes = NUMBER_OF_NODES_GEOMETRY
 593!   TOTAL_NUMBER_OF_ALL_NODES = TOTAL_NUMBER_OF_NODES + TotalNumberOfSolidNodes
 594!   MESH_NUMBER_OF_ALL_COMPONENTS = MESH_NUMBER_OF_COMPONENTS
 595! 
 596!   !Start the creation of mesh nodes
 597!   CALL CMISSNodes_Initialise(Nodes,Err)
 598!   CALL CMISSNodes_CreateStart(Region,TOTAL_NUMBER_OF_ALL_NODES,Nodes,Err)
 599!   CALL CMISSNodes_CreateFinish(Nodes,Err)
 600!   !Start the creation of the mesh
 601!   CALL CMISSMesh_Initialise(Mesh,Err)
 602!   CALL CMISSMesh_CreateStart(MeshUserNumber,Region,NUMBER_OF_DIMENSIONS,Mesh,Err)
 603!   !Set number of mesh elements
 604!   CALL CMISSMesh_NumberOfElementsSet(Mesh,TOTAL_NUMBER_OF_ELEMENTS,Err)
 605!   !Set number of mesh components
 606!   CALL CMISSMesh_NumberOfComponentsSet(Mesh,MESH_NUMBER_OF_ALL_COMPONENTS,Err)
 607!   !
 608!   CALL CMISSMeshElements_Initialise(MeshElementsGeometry,Err)
 609!   CALL CMISSMeshElements_Initialise(MeshElementsConcOne,Err)
 610!   CALL CMISSMeshElements_Initialise(MeshElementsConcTwo,Err)
 611!   CALL CMISSMeshElements_Initialise(MeshElementsConcThree,Err)
 612  MESH_COMPONENT_NUMBER_GEOMETRY=1
 613  MESH_COMPONENT_NUMBER_CONC_ONE=1
 614  MESH_COMPONENT_NUMBER_CONC_TWO=1
 615  MESH_COMPONENT_NUMBER_CONC_THREE=1
 616!   !Specify spatial mesh component
 617!   CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_GEOMETRY,BasisGeometry,MeshElementsGeometry,Err)
 618!   DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
 619!     CALL CMISSMeshElements_NodesSet(MeshElementsGeometry,ELEMENT_NUMBER,CM%M(ELEMENT_NUMBER,1:NUMBER_OF_ELEMENT_NODES_GEOMETRY),Err)
 620!   ENDDO
 621!   CALL CMISSMeshElements_CreateFinish(MeshElementsGeometry,Err)
 622!   !Specify concentration one mesh component
 623! ! !  IF(BASIS_XI_INTERPOLATION_CONC_ONE==BASIS_XI_INTERPOLATION_GEOMETRY) THEN
 624!      MeshElementsConcOne=MeshElementsGeometry
 625! ! !   ELSE
 626! ! !     MESH_COMPONENT_NUMBER_CONC_ONE=MESH_COMPONENT_NUMBER_GEOMETRY+1
 627! !      CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_CONC_ONE,BasisConcOne,MeshElementsConcOne,Err)
 628! ! !     DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
 629! ! !       CALL CMISSMeshElements_NodesSet(MeshElementsConcOne,ELEMENT_NUMBER,CM%V(ELEMENT_NUMBER, & 
 630! ! !         & 1:NUMBER_OF_ELEMENT_NODES_CONC_ONE),Err)
 631! ! !     ENDDO
 632!  !    CALL CMISSMeshElements_CreateFinish(MeshElementsConcOne,Err)
 633! ! !  ENDIF
 634! !   !Specify concentration two mesh component
 635! ! !  IF(BASIS_XI_INTERPOLATION_CONC_TWO==BASIS_XI_INTERPOLATION_GEOMETRY) THEN
 636!      MeshElementsConcTwo=MeshElementsGeometry
 637! ! !    MESH_COMPONENT_NUMBER_CONC_TWO=MESH_COMPONENT_NUMBER_GEOMETRY
 638! ! !  ELSE IF(BASIS_XI_INTERPOLATION_CONC_TWO==BASIS_XI_INTERPOLATION_CONC_ONE) THEN
 639! ! !    MeshElementsConcTwo=MeshElementsConcOne
 640! ! !    MESH_COMPONENT_NUMBER_CONC_TWO=MESH_COMPONENT_NUMBER_CONC_ONE
 641! ! !  ELSE
 642! ! !    MESH_COMPONENT_NUMBER_CONC_TWO=MESH_COMPONENT_NUMBER_CONC_ONE+1
 643! !     CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_CONC_TWO,BasisConcTwo,MeshElementsConcTwo,Err)
 644! ! !    DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
 645! ! !      CALL CMISSMeshElements_NodesSet(MeshElementsConcTwo,ELEMENT_NUMBER,CM%P(ELEMENT_NUMBER, & 
 646! ! !        & 1:NUMBER_OF_ELEMENT_NODES_CONC_TWO),Err)
 647! ! !    ENDDO
 648! !     CALL CMISSMeshElements_CreateFinish(MeshElementsConcTwo,Err)
 649! !  ! ENDIF
 650! !   !Specify concentration three mesh component
 651! ! !  IF(BASIS_XI_INTERPOLATION_CONC_THREE==BASIS_XI_INTERPOLATION_GEOMETRY) THEN
 652!      MeshElementsConcThree=MeshElementsGeometry
 653! ! !    MESH_COMPONENT_NUMBER_CONC_THREE=MESH_COMPONENT_NUMBER_GEOMETRY
 654! ! !   ELSE IF(BASIS_XI_INTERPOLATION_CONC_THREE==BASIS_XI_INTERPOLATION_CONC_ONE) THEN
 655! ! !     MeshElementsConcThree=MeshElementsConcOne
 656! ! !    MESH_COMPONENT_NUMBER_CONC_THREE=MESH_COMPONENT_NUMBER_CONC_ONE
 657! ! !  ELSE
 658! ! !    MESH_COMPONENT_NUMBER_CONC_TWO=MESH_COMPONENT_NUMBER_CONC_ONE+1
 659! !     CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_CONC_THREE,BasisConcThree,MeshElementsConcThree,Err)
 660! ! !    DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
 661! ! !      CALL CMISSMeshElements_NodesSet(MeshElementsConcTwo,ELEMENT_NUMBER,CM%P(ELEMENT_NUMBER, & 
 662! ! !        & 1:NUMBER_OF_ELEMENT_NODES_CONC_TWO),Err)
 663! ! !    ENDDO
 664! !     CALL CMISSMeshElements_CreateFinish(MeshElementsConcThree,Err)
 665! ! !  ENDIF
 666! 
 667!   !Finish the creation of the mesh
 668!   CALL CMISSMesh_CreateFinish(Mesh,Err)
 669
 670  !Start the creation of a generated mesh in the region
 671  CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
 672  CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
 673  !Set up a regular x*y*z mesh
 674  CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_REGULAR_MESH_TYPE,Err)
 675  !Set the default basis
 676  CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,BasisGeometry,Err)   
 677  !Define the mesh on the region
 678  IF(NUMBER_GLOBAL_Z_ELEMENTS==0) THEN
 679    CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/WIDTH,HEIGHT/),Err)
 680    CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS/),Err)
 681  ELSE
 682    CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/WIDTH,HEIGHT,LENGTH/),Err)
 683    CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS, &
 684      & NUMBER_GLOBAL_Z_ELEMENTS/),Err)
 685  ENDIF
 686  !Finish the creation of a generated mesh in the region
 687  CALL CMISSMesh_Initialise(Mesh,Err)
 688  CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
 689
 690  !
 691  !================================================================================================================================
 692  !
 693  !GEOMETRIC FIELD
 694  !Create a decomposition:
 695  !All mesh components share the same decomposition
 696  CALL CMISSDecomposition_Initialise(Decomposition,Err)
 697  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
 698  !Set the decomposition to be a general decomposition with the specified number of domains
 699  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
 700  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfUserDomains,Err)
 701  !Finish the decomposition
 702  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
 703
 704  !Start to create a default (geometric) field on the region
 705  CALL CMISSField_Initialise(GeometricField,Err)
 706  CALL CMISSField_CreateStart(GeometricFieldUserNumber,Region,GeometricField,Err)
 707  !Set the field type
 708  CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
 709  !Set the decomposition to use
 710  CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
 711  !Set the scaling to use
 712  CALL CMISSField_ScalingTypeSet(GeometricField,CMISS_FIELD_NO_SCALING,Err)
 713
 714  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 715    CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, & 
 716      & MESH_COMPONENT_NUMBER_GEOMETRY,Err)
 717  ENDDO
 718  !Finish creating the field
 719  CALL CMISSField_CreateFinish(GeometricField,Err)
 720
 721  !Update the geometric field parameters
 722  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricField,Err)
 723  !Update the geometric field parameters
 724!   DO NODE_NUMBER=1,NUMBER_OF_NODES_GEOMETRY
 725!     DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 726!       VALUE=CM%N(NODE_NUMBER,COMPONENT_NUMBER)
 727!       CALL CMISSField_ParameterSetUpdateNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
 728!         & CMISS_NO_GLOBAL_DERIV,NODE_NUMBER,COMPONENT_NUMBER,VALUE,Err)
 729!     ENDDO
 730!   ENDDO
 731!   CALL CMISSField_ParameterSetUpdateStart(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
 732!   CALL CMISSField_ParameterSetUpdateFinish(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
 733  !
 734  !================================================================================================================================
 735  !
 736  !EQUATIONS SETS - USING NEW ARGUMENTS TO ALLOW FOR MULTI-COMPARTMENT MODELS
 737  DO icompartment = 1,Ncompartments
 738
 739    EquationsSetFieldUserNumberDiffusion = 100_CMISSIntg+icompartment
 740    EquationsSetUserNumberDiffusion = 200_CMISSIntg+icompartment
 741    CALL CMISSField_Initialise(EquationsSetFieldDiffusion(icompartment),Err)
 742    CALL CMISSEquationsSet_Initialise(EquationsSetDiffusion(icompartment),Err)
 743    CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDiffusion,Region,GeometricField,&
 744         & CMISS_EQUATIONS_SET_CLASSICAL_FIELD_CLASS, &
 745         & CMISS_EQUATIONS_SET_DIFFUSION_EQUATION_TYPE,CMISS_EQUATIONS_SET_MULTI_COMP_TRANSPORT_DIFFUSION_SUBTYPE,&
 746         & EquationsSetFieldUserNumberDiffusion,EquationsSetFieldDiffusion(icompartment),EquationsSetDiffusion(icompartment),&
 747         & Err)
 748    !Finish creating the equations set
 749    CALL CMISSEquationsSet_CreateFinish(EquationsSetDiffusion(icompartment),Err)
 750    CALL CMISSField_ParameterSetUpdateConstant(RegionUserNumber,EquationsSetFielduserNumberDiffusion,CMISS_FIELD_U_VARIABLE_TYPE, &
 751      & CMISS_FIELD_VALUES_SET_TYPE,1,icompartment,Err)
 752    CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDiffusion(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, &
 753      & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err)
 754  END DO 
 755
 756  !-------------------------------------------------------------------------------------
 757  ! DEPENDENT FIELD: Shared
 758  CALL CMISSField_Initialise(DependentField,Err)
 759  CALL CMISSField_CreateStart(DependentFieldUserNumber,Region,DependentField,Err)
 760  CALL CMISSField_TypeSet(DependentField,CMISS_FIELD_GENERAL_TYPE,Err)  
 761  CALL CMISSField_MeshDecompositionSet(DependentField,Decomposition,Err)
 762  CALL CMISSField_GeometricFieldSet(DependentField,GeometricField,Err) 
 763  CALL CMISSField_DependentTypeSet(DependentField,CMISS_FIELD_DEPENDENT_TYPE,Err) 
 764  !Create 2N number of variables
 765  CALL CMISSField_NumberOfVariablesSet(DependentField,2*Ncompartments,Err) 
 766  !create two variables for each compartment
 767  ALLOCATE(VariableTypes(2*Ncompartments))
 768  DO num_var=1,Ncompartments
 769     VariableTypes(2*num_var-1)=CMISS_FIELD_U_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
 770     VariableTypes(2*num_var)=CMISS_FIELD_DELUDELN_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1))
 771  ENDDO
 772  CALL CMISSField_VariableTypesSet(DependentField,VariableTypes,Err) 
 773  !loop over the number of compartments
 774  DO icompartment=1,2*Ncompartments
 775    !set dimension type
 776    CALL CMISSField_DimensionSet(DependentField,VariableTypes(icompartment), &
 777       & CMISS_FIELD_SCALAR_DIMENSION_TYPE,Err)
 778    CALL CMISSField_NumberOfComponentsSet(DependentField,VariableTypes(icompartment),1,Err)
 779    CALL CMISSField_ComponentMeshComponentSet(DependentField,VariableTypes(icompartment),1, & 
 780       & MESH_COMPONENT_NUMBER_CONC_ONE,Err)
 781  ENDDO
 782  CALL CMISSField_CreateFinish(DependentField,Err)
 783  
 784  DO icompartment = 1,Ncompartments
 785    CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDiffusion(icompartment),DependentFieldUserNumber,&
 786      & DependentField,Err)
 787    CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDiffusion(icompartment),Err)
 788  ENDDO
 789
 790!   !Set the mesh component to be used by the field components.
 791! !   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 792! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDiffusionTwo,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, & 
 793! !       & MESH_COMPONENT_NUMBER_CONC_TWO,Err)
 794! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDiffusionTwo,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, & 
 795! !       & MESH_COMPONENT_NUMBER_CONC_TWO,Err)
 796! !   ENDDO
 797!   !Set the mesh component to be used by the field components.
 798! !   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 799! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDiffusionTwo,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, & 
 800! !       & MESH_COMPONENT_NUMBER_CONC_TWO,Err)
 801! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDiffusionTwo,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, & 
 802! !       & MESH_COMPONENT_NUMBER_CONC_TWO,Err)
 803! !   ENDDO
 804!   !Set the mesh component to be used by the field components.
 805! !   DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
 806! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDiffusionTwo,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, & 
 807! !       & MESH_COMPONENT_NUMBER_CONC_TWO,Err)
 808! !     CALL CMISSField_ComponentMeshComponentSet(DependentFieldDiffusionTwo,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, & 
 809! !       & MESH_COMPONENT_NUMBER_CONC_TWO,Err)
 810! !   ENDDO
 811!   !Finish the equations set dependent field variables
 812
 813  !-------------------------------------------------------------------------------------
 814  ! INITIALISE DEPENDENT FIELDS
 815!   !Initialise dependent field (concentration one components)
 816!   CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
 817!     & 1,1.0_CMISSDP,Err)
 818
 819!   !Initialise dependent field (concentration two components)
 820!   CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U2_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
 821!     & 1,INITIAL_FIELD_DIFFUSION_TWO,Err)
 822
 823!   !Initialise dependent field (concentration three components)
 824!   CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U3_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
 825!     & 1,INITIAL_FIELD_DIFFUSION_THREE,Err)
 826
 827
 828  !
 829  !================================================================================================================================
 830  !
 831  ALLOCATE(CouplingCoeffs(Ncompartments,Ncompartments))
 832  IF(Ncompartments==2)THEN
 833!     CouplingCoeffs(1,1)=1.0E-01_CMISSDP
 834! !     CouplingCoeffs(1,2)=-1.0E-04_CMISSDP
 835! !     CouplingCoeffs(2,1)=-1.0E-04_CMISSDP
 836!     CouplingCoeffs(1,2)=1.0E-01_CMISSDP
 837!     CouplingCoeffs(2,1)=1.0E-01_CMISSDP
 838!     CouplingCoeffs(2,2)=1.0E-01_CMISSDP
 839    CouplingCoeffs(1,1)=1.0E-01_CMISSDP
 840    CouplingCoeffs(1,2)=1.0E-01_CMISSDP
 841    CouplingCoeffs(2,1)=1.0E-01_CMISSDP
 842    CouplingCoeffs(2,2)=1.0E-01_CMISSDP
 843  ELSE IF(Ncompartments==3)THEN
 844    CouplingCoeffs(1,1)=1.0E-02_CMISSDP
 845    CouplingCoeffs(1,2)=1.0E-02_CMISSDP
 846    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
 847    CouplingCoeffs(2,1)=1.0E-02_CMISSDP
 848    CouplingCoeffs(2,2)=2.0E-02_CMISSDP
 849    CouplingCoeffs(2,3)=1.0E-02_CMISSDP
 850    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
 851    CouplingCoeffs(3,2)=1.0E-02_CMISSDP
 852    CouplingCoeffs(3,3)=1.0E-02_CMISSDP
 853  ELSE IF(Ncompartments==4)THEN
 854    CouplingCoeffs(1,1)=0.0E-02_CMISSDP
 855    CouplingCoeffs(1,2)=0.0E-02_CMISSDP
 856    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
 857    CouplingCoeffs(1,4)=0.0E-02_CMISSDP
 858    CouplingCoeffs(2,1)=0.0E-02_CMISSDP
 859    CouplingCoeffs(2,2)=0.0E-02_CMISSDP
 860    CouplingCoeffs(2,3)=0.0E-02_CMISSDP
 861    CouplingCoeffs(2,4)=0.0E-02_CMISSDP
 862    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
 863    CouplingCoeffs(3,2)=0.0E-02_CMISSDP
 864    CouplingCoeffs(3,3)=0.0E-02_CMISSDP
 865    CouplingCoeffs(3,4)=0.0E-02_CMISSDP
 866    CouplingCoeffs(4,1)=0.0E-02_CMISSDP
 867    CouplingCoeffs(4,2)=0.0E-02_CMISSDP
 868    CouplingCoeffs(4,3)=0.0E-02_CMISSDP
 869    CouplingCoeffs(4,4)=0.0E-02_CMISSDP
 870  ELSE IF(Ncompartments==5)THEN
 871    CouplingCoeffs(1,1)=0.0E-02_CMISSDP
 872    CouplingCoeffs(1,2)=0.0E-02_CMISSDP
 873    CouplingCoeffs(1,3)=0.0E-02_CMISSDP
 874    CouplingCoeffs(1,4)=0.0E-02_CMISSDP
 875    CouplingCoeffs(1,5)=0.0E-02_CMISSDP
 876    CouplingCoeffs(2,1)=0.0E-02_CMISSDP
 877    CouplingCoeffs(2,2)=0.0E-02_CMISSDP
 878    CouplingCoeffs(2,3)=0.0E-02_CMISSDP
 879    CouplingCoeffs(2,4)=0.0E-02_CMISSDP
 880    CouplingCoeffs(2,5)=0.0E-02_CMISSDP
 881    CouplingCoeffs(3,1)=0.0E-02_CMISSDP
 882    CouplingCoeffs(3,2)=0.0E-02_CMISSDP
 883    CouplingCoeffs(3,3)=0.0E-02_CMISSDP
 884    CouplingCoeffs(3,4)=0.0E-02_CMISSDP
 885    CouplingCoeffs(3,5)=0.0E-02_CMISSDP
 886    CouplingCoeffs(4,1)=0.0E-02_CMISSDP
 887    CouplingCoeffs(4,2)=0.0E-02_CMISSDP
 888    CouplingCoeffs(4,3)=0.0E-02_CMISSDP
 889    CouplingCoeffs(4,4)=0.0E-02_CMISSDP
 890    CouplingCoeffs(4,5)=0.0E-02_CMISSDP
 891    CouplingCoeffs(5,1)=0.0E-02_CMISSDP
 892    CouplingCoeffs(5,2)=0.0E-02_CMISSDP
 893    CouplingCoeffs(5,3)=0.0E-02_CMISSDP
 894    CouplingCoeffs(5,4)=0.0E-02_CMISSDP
 895    CouplingCoeffs(5,5)=0.0E-02_CMISSDP
 896  ELSE
 897    write(*,*) "Can't initialise coupling coefficients array."
 898  ENDIF
 899  !MATERIALS FIELDS - create the materials field
 900  !Auto-created field contains a U variable type to store the diffusion coefficient(s)
 901  !It also contains a V variable type to store the coupling coefficients 
 902  DO icompartment = 1,Ncompartments
 903    MaterialsFieldUserNumberDiffusion = 400+icompartment
 904    CALL CMISSField_Initialise(MaterialsFieldDiffusion(icompartment),Err)
 905    CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDiffusion(icompartment),MaterialsFieldUserNumberDiffusion,&
 906         & MaterialsFieldDiffusion(icompartment),Err)
 907    CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDiffusion(icompartment),Err)
 908  END DO
 909  !Initialise the coupling coefficients
 910  !Need to devise a neater way of specifying these components - e.g. specify only the upper diagonal components, and then automatically fill out the rest
 911
 912
 913  DO icompartment = 1, Ncompartments
 914    DO COMPONENT_NUMBER=1, Ncompartments
 915!       CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDiffusion(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
 916!          & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
 917        CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDiffusion(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, &
 918          & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err)
 919    END DO
 920  END DO
 921
 922!   CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDiffusionOne,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
 923!     & MaterialsFieldUserNumberDiffusionOne,POROSITY_PARAM_MAT_PROPERTIES,Err)
 924!   CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDiffusionTwo,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
 925!     & MaterialsFieldUserNumberMatPropertiesPermOverVis,PERM_OVER_VIS_PARAM_MAT_PROPERTIES,Err)
 926
 927  !
 928  !================================================================================================================================
 929  !
 930
 931  !SOURCE FIELDS
 932
 933   !create the equations set source field variables for both equations sets 
 934  DO icompartment = 1,Ncompartments
 935    SourceFieldUserNumberDiffusion=700_CMISSIntg+icompartment
 936    CALL CMISSField_Initialise(SourceFieldDiffusion(icompartment),Err)
 937    CALL CMISSEquationsSet_SourceCreateStart(EquationsSetDiffusion(icompartment),SourceFieldUserNumberDiffusion,&
 938         & SourceFieldDiffusion(icompartment),Err)
 939    CALL CMISSEquationsSet_SourceCreateFinish(EquationsSetDiffusion(icompartment),Err)
 940  END DO 
 941
 942  !Create the equations set analytic field variables
 943  DO icompartment = 1,Ncompartments
 944    AnalyticFieldUserNumberDiffusion=900_CMISSIntg+icompartment
 945    CALL CMISSField_Initialise(AnalyticFieldDiffusion(icompartment),Err)
 946
 947    CALL CMISSEquationsSet_AnalyticCreateStart(EquationsSetDiffusion(icompartment), &
 948      & CMISS_EQUATIONS_SET_MULTI_COMP_DIFFUSION_TWO_COMP_TWO_DIM,&
 949      & AnalyticFieldUserNumberDiffusion,AnalyticFieldDiffusion(icompartment),Err)
 950  
 951    !Finish the equations set analytic field variables
 952    CALL CMISSEquationsSet_AnalyticCreateFinish(EquationsSetDiffusion(icompartment),Err)
 953  END DO 
 954
 955  !
 956  !================================================================================================================================
 957  !
 958
 959  !EQUATIONS
 960  DO icompartment=1,Ncompartments
 961    CALL CMISSEquations_Initialise(EquationsDiffusion(icompartment),Err)
 962    CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDiffusion(icompartment),EquationsDiffusion(icompartment),Err)
 963    CALL CMISSEquations_SparsityTypeSet(EquationsDiffusion(icompartment),CMISS_EQUATIONS_SPARSE_MATRICES,Err)
 964    CALL CMISSEquations_OutputTypeSet(EquationsDiffusion(icompartment),EQUATIONS_DIFFUSION_OUTPUT,Err)
 965    CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetDiffusion(icompartment),Err)
 966  ENDDO
 967  !
 968  !================================================================================================================================
 969  !
 970  !PROBLEMS
 971
 972  !Start the creation of a problem.
 973  CALL CMISSProblem_Initialise(Problem,Err)
 974  CALL CMISSControlLoop_Initialise(ControlLoop,Err)
 975  CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
 976  !Set the problem to be a coupled diffusion-diffusion problem
 977  CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_MULTI_PHYSICS_CLASS,CMISS_PROBLEM_MULTI_COMPARTMENT_TRANSPORT_TYPE, &
 978    & CMISS_PROBLEM_STANDARD_MULTI_COMPARTMENT_TRANSPORT_SUBTYPE,Err)
 979  !Finish the creation of a problem.
 980  CALL CMISSProblem_CreateFinish(Problem,Err)
 981  !Start the creation of the problem control loop
 982  CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
 983  !Get the control loop
 984  CALL CMISSProblem_ControlLoopGet(Problem,CMISS_CONTROL_LOOP_NODE,ControlLoop,Err)
 985  !Set the times
 986  CALL CMISSControlLoop_TimesSet(ControlLoop,LINEAR_SOLVER_DIFFUSION_START_TIME,&
 987    & LINEAR_SOLVER_DIFFUSION_STOP_TIME,LINEAR_SOLVER_DIFFUSION_TIME_INCREMENT,Err)
 988  !Set the output timing
 989  CALL CMISSControlLoop_TimeOutputSet(ControlLoop,LINEAR_SOLVER_DIFFUSION_OUTPUT_FREQUENCY,Err)
 990  !Finish creating the problem control loop
 991  CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
 992
 993  !
 994  !================================================================================================================================
 995  !
 996  !SOLVERS
 997
 998  !Start the creation of the problem solvers
 999  CALL CMISSSolver_Initialise(SolverDiffusion,Err)
1000  CALL CMISSSolver_Initialise(LinearSolverDiffusion,Err)
1001  CALL CMISSProblem_SolversCreateStart(Problem,Err)
1002
1003  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,SolverDiffusionUserNumber,SolverDiffusion,Err)
1004  !Set the output type
1005  CALL CMISSSolver_OutputTypeSet(SolverDiffusion,LINEAR_SOLVER_DIFFUSION_OUTPUT_TYPE,Err)
1006  !Set the solver settings
1007!  IF(LINEAR_SOLVER_DIFFUSION_ONE_DIRECT_FLAG) THEN
1008!    CALL CMISSSolver_LinearTypeSet(LinearSolverDiffusionOne,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
1009!    CALL CMISSSolver_LibraryTypeSet(LinearSolverDiffusionOne,CMISS_SOLVER_MUMPS_LIBRARY,Err)
1010!  ELSE
1011!    CALL CMISSSolver_LinearTypeSet(LinearSolverDiffusionOne,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
1012   CALL CMISSSolver_DynamicLinearSolverGet(SolverDiffusion,LinearSolverDiffusion,Err)
1013    CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverDiffusion,MAXIMUM_ITERATIONS,Err)
1014!    CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverDiffusionOne,DIVERGENCE_TOLERANCE,Err)
1015!    CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverDiffusionOne,RELATIVE_TOLERANCE,Err)
1016!    CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverDiffusionOne,ABSOLUTE_TOLERANCE,Err)
1017!    CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverDiffusionOne,RESTART_VALUE,Err)
1018!  ENDIF
1019  !Finish the creation of the problem solver
1020  CALL CMISSProblem_SolversCreateFinish(Problem,Err)
1021
1022  !
1023  !================================================================================================================================
1024  !
1025  !SOLVER EQUATIONS
1026  !Start the creation of the problem solver equations
1027  CALL CMISSSolver_Initialise(SolverDiffusion,Err)
1028  CALL CMISSSolverEquations_Initialise(SolverEquationsDiffusion,Err)
1029  CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
1030  !
1031  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,SolverDiffusionUserNumber,SolverDiffusion,Err)
1032  CALL CMISSSolver_SolverEquationsGet(SolverDiffusion,SolverEquationsDiffusion,Err)
1033  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsDiffusion,CMISS_SOLVER_SPARSE_MATRICES,Err)
1034
1035  DO icompartment=1,Ncompartments
1036    CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsDiffusion,EquationsSetDiffusion(icompartment),&
1037         & EquationsSetIndex,Err)
1038  ENDDO
1039
1040  !Finish the creation of the problem solver equations
1041  CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
1042
1043  !
1044  !================================================================================================================================
1045  !
1046
1047  !BOUNDARY CONDITIONS
1048  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsDiffusion,BoundaryConditionsDiffusion,Err)
1049  CALL CMISSSolverEquations_BoundaryConditionsAnalytic(SolverEquationsDiffusion,Err)
1050  CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsDiffusion,Err)
1051
1052  !
1053  !================================================================================================================================
1054  !
1055  !RUN SOLVERS
1056
1057  !Turn off PETSc error handling
1058  !CALL PETSC_ERRORHANDLING_SET_ON(ERR,ERROR,*999)
1059
1060  !Solve the problem
1061  CALL CMISSProblem_Solve(Problem,Err)
1062  !
1063  !================================================================================================================================
1064  !
1065  !OUTPUT
1066
1067  Call CMISSAnalyticAnalysisOutput(DependentField,"MultiCompDiffusionAnalytics_2Comp_2D_x40_y40_L_T1",Err)
1068
1069
1070  EXPORT_FIELD_IO=.TRUE.
1071  IF(EXPORT_FIELD_IO) THEN
1072    WRITE(*,'(A)') "Exporting fields..."
1073!     CALL CMISSFields_Initialise(Fields,Err)
1074!     CALL CMISSFields_Create(Region,Fields,Err)
1075!     CALL CMISSFields_NodesExport(Fields,"MonolithicMultiCompDiffusionTest","FORTRAN",Err)
1076!     CALL CMISSFields_ElementsExport(Fields,"MonolithicMultiCompDiffusionTest","FORTRAN",Err)
1077!     CALL CMISSFields_Finalise(Fields,Err)
1078    WRITE(*,'(A)') "Field exported!"
1079    CALL FieldmlOutput_InitializeInfo( Region, Mesh, 3, outputDirectory, basename, fieldmlInfo, err )
1080
1081    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".geometric", region, mesh, GeometricField, & 
1082      & CMISS_FIELD_U_VARIABLE_TYPE, err )
1083
1084!     CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".dependent", region, mesh, DependentField, err )
1085
1086    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".dependent_U", region, mesh, DependentField, &
1087      & CMISS_FIELD_U_VARIABLE_TYPE, err )
1088
1089    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".dependent_V", region, mesh, DependentField, &
1090      & CMISS_FIELD_V_VARIABLE_TYPE, err )
1091    
1092    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".equations_set_field_1", region, mesh, EquationsSetFieldDiffusion(1), &
1093      & CMISS_FIELD_U_VARIABLE_TYPE, err )
1094
1095    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".equations_set_field_2", region, mesh, EquationsSetFieldDiffusion(2), &
1096      & CMISS_FIELD_U_VARIABLE_TYPE, err )
1097
1098!     CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".equations_set_field_3", region, mesh, EquationsSetFieldDiffusion(3), &
1099!       & CMISS_FIELD_U_VARIABLE_TYPE, err )
1100
1101    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".source_1", region,mesh,SourceFieldDiffusion(1), &
1102      & CMISS_FIELD_U_VARIABLE_TYPE,&
1103      & err )
1104
1105    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".source_2", region,mesh,SourceFieldDiffusion(2), &
1106      & CMISS_FIELD_U_VARIABLE_TYPE,&
1107      & err )
1108
1109    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".materials_1", region, mesh, MaterialsFieldDiffusion(1), &
1110      & CMISS_FIELD_U_VARIABLE_TYPE,err )
1111
1112    CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".materials_2", region, mesh, MaterialsFieldDiffusion(2), &
1113      & CMISS_FIELD_U_VARIABLE_TYPE,err )
1114
1115    !CALL FieldmlOutput_AddField( fieldmlInfo, baseName//".analytic", region, mesh, AnalyticField, err )
1116    
1117    CALL FieldmlOutput_Write( fieldmlInfo, outputFilename, err )
1118    
1119    CALL FieldmlUtil_FinalizeInfo( fieldmlInfo )
1120  ENDIF
1121
1122  !Finialise CMISS
1123!   CALL CMISSFinalise(Err)
1124
1125  WRITE(*,'(A)') "Program successfully completed."
1126  
1127  STOP
1128
1129END PROGRAM MONOLITHICSCHEMETESTFIELDMLEXAMPLE