PageRenderTime 89ms CodeModel.GetById 16ms app.highlight 64ms RepoModel.GetById 1ms app.codeStats 0ms

/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/IncompressibleElasticityDrivenDarcyIO/src/IncompressibleElasticityDrivenDarcyIOExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 1137 lines | 571 code | 184 blank | 382 comment | 0 complexity | fa0189050314fecf43370947ff8449ab MD5 | raw file

Large files files are truncated, but you can click here to view the full 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 FINITEELASTICITYDARCYIOEXAMPLE
 57
 58  !
 59  !================================================================================================================================
 60  !
 61
 62  !PROGRAM LIBRARIES
 63
 64  USE OPENCMISS
 65!   USE FLUID_MECHANICS_IO_ROUTINES
 66  USE IOSTUFF
 67  USE MPI
 68
 69#ifdef WIN32
 70  USE IFQWINCMISS
 71#endif
 72
 73  !
 74  !================================================================================================================================
 75  !
 76
 77  !PROGRAM VARIABLES AND TYPES
 78
 79  IMPLICIT NONE
 80
 81  !Test program parameters
 82
 83  REAL(CMISSDP), PARAMETER :: Y_DIM=1.0_CMISSDP
 84  REAL(CMISSDP), PARAMETER :: X_DIM=1.0_CMISSDP
 85  REAL(CMISSDP), PARAMETER :: Z_DIM=1.0_CMISSDP
 86
 87  INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=1
 88  INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=2
 89  INTEGER(CMISSIntg), PARAMETER :: CubicBasisUserNumber=3
 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 :: MaterialsFieldUserNumberDarcy=8
 97  INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberDarcy=12
 98  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=14
 99  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumberDarcy=22
100
101  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1
102  INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2
103  INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1
104  INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1
105  INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1
106  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPorosity=1
107  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPermOverVis=2
108
109  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1
110  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
111
112  !Program types
113
114  !Program variables
115
116  INTEGER(CMISSIntg) :: NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS,NUMBER_GLOBAL_Z_ELEMENTS
117
118!   INTEGER(CMISSIntg) :: MPI_IERROR
119  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,NumberOfDomains,ComputationalNodeNumber
120
121  INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
122
123  INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
124  INTEGER(CMISSIntg) :: RESTART_VALUE
125
126  INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT
127  INTEGER(CMISSIntg) :: COMPONENT_NUMBER
128
129  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY
130  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE
131  INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE
132  INTEGER(CMISSIntg) :: LINEAR_SOLVER_MAT_PROPERTIES_OUTPUT_TYPE
133
134  REAL(CMISSDP) :: GEOMETRY_TOLERANCE
135  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID
136  REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4)
137  REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
138  REAL(CMISSDP) :: RELATIVE_TOLERANCE
139  REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
140  REAL(CMISSDP) :: LINESEARCH_ALPHA
141  REAL(CMISSDP) :: VALUE
142  REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY
143
144  LOGICAL :: EXPORT_FIELD_IO
145  LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG
146
147  !CMISS variables
148
149  !Regions
150  TYPE(CMISSRegionType) :: Region
151  TYPE(CMISSRegionType) :: WorldRegion
152  !Coordinate systems
153  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem
154  TYPE(CMISSCoordinateSystemType) :: WorldCoordinateSystem
155  !Basis
156!   TYPE(CMISSBasisType) :: CubicBasis, QuadraticBasis, LinearBasis, Bases(2)
157  TYPE(CMISSBasisType),allocatable :: Bases(:)
158  !Meshes
159  TYPE(CMISSMeshType) :: Mesh
160  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
161
162  !Decompositions
163  TYPE(CMISSDecompositionType) :: Decomposition
164  !Fields
165  TYPE(CMISSFieldsType) :: Fields
166  !Field types
167  TYPE(CMISSFieldType) :: GeometricField
168  TYPE(CMISSFieldType) :: MaterialsFieldDarcy
169  TYPE(CMISSFieldType) :: EquationsSetFieldDarcy
170  !Boundary conditions
171  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy
172  !Equations sets
173  TYPE(CMISSEquationsSetType) :: EquationsSetDarcy
174  !Equations
175  TYPE(CMISSEquationsType) :: EquationsDarcy
176  !Problems
177  TYPE(CMISSProblemType) :: Problem
178  !Control loops
179  TYPE(CMISSControlLoopType) :: ControlLoop
180  !Solvers
181  TYPE(CMISSSolverType) :: DynamicSolverDarcy
182  TYPE(CMISSSolverType) :: LinearSolverDarcy
183!   TYPE(CMISSSolverType) :: LinearSolverSolid
184  !Solver equations
185  TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy
186  ! nodes and elements
187  TYPE(CMISSNodesType) :: Nodes
188  TYPE(CMISSMeshElementsType),allocatable :: Elements(:)
189
190  !Other variables
191  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face1Nodes(:),Face2Nodes(:)
192  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face3Nodes(:),Face4Nodes(:)
193  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face5Nodes(:),Face6Nodes(:)
194  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face7Nodes(:),Face8Nodes(:)
195  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face9Nodes(:),Face10Nodes(:)
196  INTEGER(CMISSIntg),ALLOCATABLE,TARGET :: Face11Nodes(:),Face12Nodes(:)
197  INTEGER(CMISSIntg) :: FaceXi(6)
198  INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
199  REAL(CMISSDP) :: XCoord,YCoord,ZCoord
200  LOGICAL :: X_FIXED,Y_FIXED !,X_OKAY,Y_OKAY
201
202#ifdef WIN32
203  !Quickwin type
204  LOGICAL :: QUICKWIN_STATUS=.FALSE.
205  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
206#endif
207
208  !Generic CMISS variables
209
210  INTEGER(CMISSIntg) :: EquationsSetIndex
211  INTEGER(CMISSIntg) :: Err
212
213
214  INTEGER(CMISSIntg) :: DIAG_LEVEL_LIST(5)
215!   CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(8) !,TIMING_ROUTINE_LIST(1)
216  CHARACTER(LEN=255) :: DIAG_ROUTINE_LIST(1) !,TIMING_ROUTINE_LIST(1)
217
218  !
219  !--------------------------------------------------------------------------------------------------------------------------------
220  !
221
222  !Program variables and types (finite elasticity part)
223
224  !Test program parameters
225
226  INTEGER(CMISSIntg) :: SolidMeshComponenetNumber
227
228  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidUserNumber=1
229  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfVariables=1
230  INTEGER(CMISSIntg), PARAMETER :: FieldGeometrySolidNumberOfComponents=3
231
232  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidUserNumber=2
233  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfVariables=1
234  INTEGER(CMISSIntg), PARAMETER :: FieldFibreSolidNumberOfComponents=3
235
236  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidUserNumber=3
237  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfVariables=1
238  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialSolidNumberOfComponents=3
239
240  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidUserNumber=4
241  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfVariables=4
242  INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4
243  INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4  !(u,v,w,m)
244
245  INTEGER(CMISSIntg), PARAMETER :: EquationSetSolidUserNumber=1
246  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldSolidUserNumber=25
247
248  INTEGER(CMISSIntg), PARAMETER :: SolidDisplMeshComponentNumber=1
249  INTEGER(CMISSIntg), PARAMETER :: SolidLagrMultMeshComponentNumber=2
250  INTEGER(CMISSIntg), PARAMETER :: SolidGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
251
252  INTEGER(CMISSIntg), PARAMETER :: DarcyVelMeshComponentNumber=SolidLagrMultMeshComponentNumber
253  INTEGER(CMISSIntg), PARAMETER :: DarcyMassIncreaseMeshComponentNumber=SolidLagrMultMeshComponentNumber
254!   INTEGER(CMISSIntg), PARAMETER :: DarcyGeometryMeshComponentNumber=SolidDisplMeshComponentNumber
255
256  INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=32
257  !Program types
258  !Program variables
259
260  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME
261  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME
262  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA
263  REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
264
265  !CMISS variables
266
267  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsSolid
268  TYPE(CMISSEquationsType) :: EquationsSolid
269  TYPE(CMISSEquationsSetType) :: EquationsSetSolid
270  TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid
271  TYPE(CMISSFieldType) :: DependentFieldSolid,EquationsSetFieldSolid
272  TYPE(CMISSSolverType) :: SolverSolid
273  TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid
274
275  !End - Program variables and types (finite elasticity part)
276
277  !
278  !--------------------------------------------------------------------------------------------------------------------------------
279  !
280
281
282#ifdef WIN32
283  !Initialise QuickWin
284  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
285  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
286  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
287  !Set the window parameters
288  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
289  !If attempt fails set with system estimated values
290  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
291#endif
292
293  !
294  !================================================================================================================================
295  !
296  NUMBER_GLOBAL_X_ELEMENTS=1
297  NUMBER_GLOBAL_Y_ELEMENTS=1
298  NUMBER_GLOBAL_Z_ELEMENTS=3
299
300  IF(NUMBER_GLOBAL_Z_ELEMENTS==0)THEN
301    NUMBER_OF_DIMENSIONS=2
302  ELSE
303    NUMBER_OF_DIMENSIONS=3
304  ENDIF
305  !PROBLEM CONTROL PANEL
306
307!   BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION
308  BASIS_XI_INTERPOLATION_SOLID=CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION
309  !Set geometric tolerance
310  GEOMETRY_TOLERANCE = 1.0E-12_CMISSDP
311  !Set initial values
312  INITIAL_FIELD_DARCY(1)=0.0_CMISSDP
313  INITIAL_FIELD_DARCY(2)=0.0_CMISSDP
314  INITIAL_FIELD_DARCY(3)=0.0_CMISSDP
315  INITIAL_FIELD_DARCY(4)=0.0_CMISSDP
316  !Set material parameters
317  POROSITY_PARAM_DARCY=0.1_CMISSDP
318  PERM_OVER_VIS_PARAM_DARCY=1.0_CMISSDP
319  !Set output parameter
320  !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
321  DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT
322  LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_SOLVER_OUTPUT
323  !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
324  EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
325
326  !Set time parameter
327  DYNAMIC_SOLVER_DARCY_START_TIME=0.0_CMISSDP
328  DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP
329  DYNAMIC_SOLVER_DARCY_STOP_TIME=2_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT
330  DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP
331  !Set result output parameter
332  DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1
333  !Set solver parameters
334  LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE.
335  RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
336  ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
337  DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5
338  MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000
339  RESTART_VALUE=30_CMISSIntg !default: 30
340  LINESEARCH_ALPHA=1.0_CMISSDP
341
342
343  !
344  !================================================================================================================================
345  !
346
347  !INITIALISE OPENCMISS
348
349  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
350
351  CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
352
353  !
354  !================================================================================================================================
355  !
356
357  !Set diagnostics
358
359  DIAG_LEVEL_LIST(1)=1
360  DIAG_LEVEL_LIST(2)=2
361  DIAG_LEVEL_LIST(3)=3
362  DIAG_LEVEL_LIST(4)=4
363  DIAG_LEVEL_LIST(5)=5
364
365  DIAG_ROUTINE_LIST(1)="WRITE_IP_INFO"
366!   DIAG_ROUTINE_LIST(2)="FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR"
367
368  !CMISS_ALL_DIAG_TYPE/CMISS_IN_DIAG_TYPE/CMISS_FROM_DIAG_TYPE
369  CALL CMISSDiagnosticsSetOn(CMISS_IN_DIAG_TYPE,DIAG_LEVEL_LIST,"Diagnostics",DIAG_ROUTINE_LIST,Err)
370
371  !CMISS_ALL_TIMING_TYPE/CMISS_IN_TIMING_TYPE/CMISS_FROM_TIMING_TYPE
372  !TIMING_ROUTINE_LIST(1)="PROBLEM_FINITE_ELEMENT_CALCULATE"
373  !CALL TIMING_SET_ON(IN_TIMING_TYPE,.TRUE.,"",TIMING_ROUTINE_LIST,ERR,ERROR,*999)
374
375  !
376  !================================================================================================================================
377  !
378
379  !Get the number of computational nodes and this computational node number
380  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
381  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
382
383  NumberOfDomains = NumberOfComputationalNodes
384  write(*,*) "NumberOfDomains = ",NumberOfDomains
385
386  !
387  !================================================================================================================================
388  !
389
390  !COORDINATE SYSTEM
391
392  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
393  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
394  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NUMBER_OF_DIMENSIONS,Err)
395  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
396
397  !
398  !================================================================================================================================
399  !
400
401  !REGION
402  !For a volume-coupled problem, solid and fluid are based in the same region
403
404  CALL CMISSRegion_Initialise(Region,Err)
405  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
406  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
407  CALL CMISSRegion_CreateFinish(Region,Err)
408
409  !
410  !================================================================================================================================
411  !
412
413call READ_MESH('input/example-mesh-3els',MeshUserNumber,Region, Mesh,Bases,Nodes,Elements)
414!   !BASES
415!   !Define basis functions
416!   CALL CMISSBasis_Initialise(LinearBasis,Err)
417!   CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
418!   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
419!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
420!   !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err)
421!   CALL CMISSBasis_CreateFinish(LinearBasis,Err)
422! 
423!   CALL CMISSBasis_Initialise(QuadraticBasis,Err)
424!   CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
425!   CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
426!     & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
427!   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
428!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
429!   !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err)
430!   CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
431! 
432!   CALL CMISSBasis_Initialise(CubicBasis,Err)
433!   CALL CMISSBasis_CreateStart(CubicBasisUserNumber,CubicBasis,Err)
434!   CALL CMISSBasis_InterpolationXiSet(CubicBasis,(/CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION, &
435!     & CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_CUBIC_LAGRANGE_INTERPOLATION/),Err)
436!   CALL CMISSBasis_QuadratureNumberOfGaussXiSet(CubicBasis, &
437!     & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err)
438!   !CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(CubicBasis,.true.,Err) !Enable 3D interpolation on faces
439!   CALL CMISSBasis_CreateFinish(CubicBasis,Err)
440! 
441!   !LinearBasis/QuadraticBasis/CubicBasis
442!   Bases(1)=QuadraticBasis
443!   Bases(2)=LinearBasis
444! 
445!   !Start the creation of a generated mesh in the region
446!   CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
447!   CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
448!   CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_REGULAR_MESH_TYPE,Err)
449!   CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,Bases,Err)
450!   IF(NUMBER_GLOBAL_Z_ELEMENTS==0) THEN
451!     CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/X_DIM,Y_DIM/),Err)
452!     CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS/),Err)
453!   ELSE
454!     CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/X_DIM,Y_DIM,Z_DIM/),Err)
455!     CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NUMBER_GLOBAL_X_ELEMENTS,NUMBER_GLOBAL_Y_ELEMENTS, &
456!       & NUMBER_GLOBAL_Z_ELEMENTS/),Err)
457!   ENDIF
458!   CALL CMISSMesh_Initialise(Mesh,Err)
459!   CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
460! -------------------------------------------------------------------------------
461  
462
463  !GEOMETRIC FIELD
464
465  !Create a decomposition:
466  CALL CMISSDecomposition_Initialise(Decomposition,Err)
467  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
468  !Set the decomposition to be a general decomposition with the specified number of domains
469  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
470  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
471  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
472
473  CALL CMISSField_Initialise(GeometricField,Err)
474  CALL CMISSField_CreateStart(GeometricFieldUserNumber,Region,GeometricField,Err)
475  CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
476  CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
477  CALL CMISSField_NumberOfVariablesSet(GeometricField,FieldGeometryNumberOfVariables,Err)
478  CALL CMISSField_NumberOfComponentsSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,Err)  
479  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
480  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
481  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
482  CALL CMISSField_CreateFinish(GeometricField,Err)
483CALL READ_NODES('input/example-nodes-3els',GeometricField)
484!   CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricField,Err)
485
486  !--------------------------------------------------------------------------------------------------------------------------------
487  ! Solid
488
489  !Create a decomposition
490
491  !Create a field to put the geometry (defualt is geometry)
492
493  SolidMeshComponenetNumber = SolidGeometryMeshComponentNumber
494
495  CALL CMISSField_Initialise(GeometricFieldSolid,Err)
496  CALL CMISSField_CreateStart(FieldGeometrySolidUserNumber,Region,GeometricFieldSolid,Err)
497  CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err)
498  CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err)
499  CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometrySolidNumberOfVariables,Err)
500  CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometrySolidNumberOfComponents,Err)
501  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidMeshComponenetNumber,Err)
502  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidMeshComponenetNumber,Err)
503  CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidMeshComponenetNumber,Err)
504  CALL CMISSField_CreateFinish(GeometricFieldSolid,Err)
505  !Set the mesh component to be used by the field components.
506CALL READ_NODES('input/example-nodes-3els',GeometricFieldSolid)
507!  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err)
508
509  !Create a fibre field and attach it to the geometric field
510  CALL CMISSField_Initialise(FibreFieldSolid,Err)
511  CALL CMISSField_CreateStart(FieldFibreSolidUserNumber,Region,FibreFieldSolid,Err)
512  CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err)
513  CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err)
514  CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err)
515  CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreSolidNumberOfVariables,Err)
516  CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreSolidNumberOfComponents,Err)
517  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
518  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
519  CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
520  CALL CMISSField_CreateFinish(FibreFieldSolid,Err)
521
522  ! end Solid
523  !--------------------------------------------------------------------------------------------------------------------------------
524
525  !
526  !================================================================================================================================
527  !
528
529  !EQUATIONS SETS
530
531  !Create the equations set for ALE Darcy
532  CALL CMISSField_Initialise(EquationsSetFieldDarcy,Err)
533  CALL CMISSEquationsSet_Initialise(EquationsSetDarcy,Err)
534  CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricField,CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
535    & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,&
536    & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy,EquationsSetDarcy,Err)
537  CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy,Err)
538
539  !Create the equations set for the solid
540  CALL CMISSField_Initialise(EquationsSetFieldSolid,Err)
541  CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err)
542  CALL CMISSEquationsSet_CreateStart(EquationSetSolidUserNumber,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
543    & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,&
544    & EquationsSetFieldSolidUserNumber,EquationsSetFieldSolid,EquationsSetSolid,Err)
545  CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err)
546
547  !--------------------------------------------------------------------------------------------------------------------------------
548  ! Solid Materials Field
549
550  !Create a material field and attach it to the geometric field
551  CALL CMISSField_Initialise(MaterialFieldSolid,Err)
552  !
553  CALL CMISSField_CreateStart(FieldMaterialSolidUserNumber,Region,MaterialFieldSolid,Err)
554  !
555  CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err)
556  CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err)
557  CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err)
558  CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialSolidNumberOfVariables,Err)
559  CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialSolidNumberOfComponents,Err)
560  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidGeometryMeshComponentNumber,Err)
561  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidGeometryMeshComponentNumber,Err)
562  CALL CMISSField_ComponentMeshComponentSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidGeometryMeshComponentNumber,Err)
563  !
564  CALL CMISSField_CreateFinish(MaterialFieldSolid,Err)
565
566  !Set material parameters
567  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
568    & 2.0_CMISSDP,Err)
569!   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,2.0e3_CMISSDP,Err)
570  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2, &
571    & 6.0_CMISSDP,Err)
572!   CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,33.0_CMISSDP,Err)
573  CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3, &
574    & 10.0_CMISSDP,Err)
575
576
577  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialSolidUserNumber,MaterialFieldSolid,Err)
578  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err)
579
580  ! end Solid
581  !--------------------------------------------------------------------------------------------------------------------------------
582
583
584  !
585  !================================================================================================================================
586  !
587
588  !DEPENDENT FIELDS
589
590  !--------------------------------------------------------------------------------------------------------------------------------
591  ! Solid
592
593  !Create a dependent field with four variables (U, DelUDelN = solid, V, DelVDelN = Darcy) and four components
594  CALL CMISSField_Initialise(DependentFieldSolid,Err)
595  !
596  CALL CMISSField_CreateStart(FieldDependentSolidUserNumber,Region,DependentFieldSolid,Err)
597  !
598  CALL CMISSField_TypeSet(DependentFieldSolid,CMISS_FIELD_GENERAL_TYPE,Err)
599  CALL CMISSField_MeshDecompositionSet(DependentFieldSolid,Decomposition,Err)
600  CALL CMISSField_GeometricFieldSet(DependentFieldSolid,GeometricFieldSolid,Err)
601  CALL CMISSField_DependentTypeSet(DependentFieldSolid,CMISS_FIELD_DEPENDENT_TYPE,Err)
602  CALL CMISSField_NumberOfVariablesSet(DependentFieldSolid,FieldDependentSolidNumberOfVariables,Err)
603  CALL CMISSField_VariableTypesSet(DependentFieldSolid,(/CMISS_FIELD_U_VARIABLE_TYPE, &
604    & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err)
605  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err)
606  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
607    & FieldDependentSolidNumberOfComponents,Err)
608  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err)
609  CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE, &
610    & FieldDependentFluidNumberOfComponents,Err)
611  !
612  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidDisplMeshComponentNumber,Err)
613  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidDisplMeshComponentNumber,Err)
614  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidDisplMeshComponentNumber,Err)
615  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4, &
616    & CMISS_FIELD_NODE_BASED_INTERPOLATION, &
617    & Err)
618!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
619!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
620  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidLagrMultMeshComponentNumber,Err)
621  !
622  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, &
623    & SolidDisplMeshComponentNumber, &
624    & Err)
625  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, &
626    & SolidDisplMeshComponentNumber, &
627    & Err)
628  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, &
629    & SolidDisplMeshComponentNumber, &
630    & Err)
631  CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
632    & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err)
633!   CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
634!     & CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err)
635!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err)
636  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, &
637    & SolidLagrMultMeshComponentNumber, &
638    & Err)
639
640  !For this equation type, MESH_COMPONENT_NUMBER_PRESSURE is actually the mass increase component as the pressure is taken from the solid equations
641  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber,Err)
642  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber,Err)
643  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber,Err)
644!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
645  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4, &
646    & DarcyMassIncreaseMeshComponentNumber, &
647    & Err)
648  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber, &
649    & Err)
650  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber, &
651    & Err)
652  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber, &
653    & Err)
654!   CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err)
655  CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4, &
656    & DarcyMassIncreaseMeshComponentNumber,Err)
657
658  CALL CMISSField_CreateFinish(DependentFieldSolid,Err)
659  !
660  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentSolidUserNumber,DependentFieldSolid,Err)
661  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err)
662
663  ! end Solid
664  !--------------------------------------------------------------------------------------------------------------------------------
665
666
667  !Create the equations set dependent field variables for ALE Darcy
668  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,FieldDependentSolidUserNumber, & ! ??? UserNumber ???
669    & DependentFieldSolid,Err)
670  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy,Err)
671
672  !Initialise dependent field (velocity components,mass increase)
673  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS+1
674    CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
675      & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err)
676  ENDDO
677
678
679  !
680  !================================================================================================================================
681  !
682
683call READ_FIELD('input/example-field-darcy',MaterialsFieldUserNumberDarcy,Region,GeometricField,MaterialsFieldDarcy)
684CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy,MaterialsFieldUserNumberDarcy,MaterialsFieldDarcy,Err)
685CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy,Err)
686!   !MATERIALS FIELDS
687! 
688!   !Create the equations set materials field variables for ALE Darcy
689!   CALL CMISSField_Initialise(MaterialsFieldDarcy,Err)
690!   CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy,MaterialsFieldUserNumberDarcy, &
691!     & MaterialsFieldDarcy,Err)
692!   !Finish the equations set materials field variables
693!   CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy,Err)
694!   CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
695!     & MaterialsFieldUserNumberDarcyPorosity,POROSITY_PARAM_DARCY,Err)
696!   CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
697!     & MaterialsFieldUserNumberDarcyPermOverVis,PERM_OVER_VIS_PARAM_DARCY,Err)
698
699  !
700  !================================================================================================================================
701  !
702
703  !EQUATIONS SET EQUATIONS
704
705  !Darcy
706  CALL CMISSEquations_Initialise(EquationsDarcy,Err)
707  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy,EquationsDarcy,Err)
708  CALL CMISSEquations_SparsityTypeSet(EquationsDarcy,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
709  CALL CMISSEquations_OutputTypeSet(EquationsDarcy,EQUATIONS_DARCY_OUTPUT,Err)
710  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetDarcy,Err)
711
712  !Solid
713  CALL CMISSEquations_Initialise(EquationsSolid,Err)
714  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetSolid,EquationsSolid,Err)
715  CALL CMISSEquations_SparsityTypeSet(EquationsSolid,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
716  CALL CMISSEquations_OutputTypeSet(EquationsSolid,CMISS_EQUATIONS_NO_OUTPUT,Err)
717  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetSolid,Err)
718
719  !
720  !================================================================================================================================
721  !
722
723  !--------------------------------------------------------------------------------------------------------------------------------
724  ! Solid
725
726  !Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure
727  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
728    & CMISS_FIELD_VALUES_SET_TYPE, &
729    & 1,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,Err)
730  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
731    & CMISS_FIELD_VALUES_SET_TYPE, &
732    & 2,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,Err)
733  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, &
734    & CMISS_FIELD_VALUES_SET_TYPE, &
735    & 3,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,Err)
736  CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4, &
737    & 0.0_CMISSDP, &
738    & Err)
739
740
741  !
742  !================================================================================================================================
743  !
744
745  !PROBLEMS
746
747  CALL CMISSProblem_Initialise(Problem,Err)
748  CALL CMISSControlLoop_Initialise(ControlLoop,Err)
749  CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
750  CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_MULTI_PHYSICS_CLASS,CMISS_PROBLEM_FINITE_ELASTICITY_DARCY_TYPE, &
751    & CMISS_PROBLEM_QUASISTATIC_ELASTICITY_TRANSIENT_DARCY_SUBTYPE,Err)
752  CALL CMISSProblem_CreateFinish(Problem,Err)
753
754  CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
755  CALL CMISSProblem_ControlLoopGet(Problem,CMISS_CONTROL_LOOP_NODE,ControlLoop,Err)
756!   CALL CMISSControlLoop_MaximumIterationsSet(ControlLoop,2,Err)
757  CALL CMISSControlLoop_TimesSet(ControlLoop,DYNAMIC_SOLVER_DARCY_START_TIME,DYNAMIC_SOLVER_DARCY_STOP_TIME, &
758    & DYNAMIC_SOLVER_DARCY_TIME_INCREMENT,Err)
759  CALL CMISSControlLoop_TimeOutputSet(ControlLoop,DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY,Err)
760!   CALL CMISSControlLoop_OutputTypeSet(ControlLoop,CMISS_CONTROL_LOOP_PROGRESS_OUTPUT,Err)
761  CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
762
763  !
764  !================================================================================================================================
765  !
766
767  !SOLVERS
768
769  CALL CMISSSolver_Initialise(SolverSolid,Err)
770  CALL CMISSSolver_Initialise(DynamicSolverDarcy,Err)
771  CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
772
773  CALL CMISSProblem_SolversCreateStart(Problem,Err)
774
775  ! Solid
776  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
777    & SolverSolidIndex,SolverSolid,Err)
778  CALL CMISSSolver_OutputTypeSet(SolverSolid,CMISS_SOLVER_PROGRESS_OUTPUT,Err)
779!   CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_FD_CALCULATED,Err)
780  CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_EQUATIONS_CALCULATED,Err)
781
782  CALL CMISSSolver_NewtonAbsoluteToleranceSet(SolverSolid,ABSOLUTE_TOLERANCE,Err)
783  CALL CMISSSolver_NewtonRelativeToleranceSet(SolverSolid,RELATIVE_TOLERANCE,Err)
784  CALL CMISSSolver_NewtonMaximumIterationsSet(SolverSolid,MAXIMUM_ITERATIONS,Err)
785
786!   CALL CMISSSolverNonLinearTypeSet(SolverSolid,CMISS_SOLVER_NONLINEAR_NEWTON,Err)
787!   CALL CMISSSolver_LibraryTypeSet(SolverSolid,CMISS_SOLVER_PETSC_LIBRARY,Err)
788
789!   CALL CMISSSolver_NewtonLinearSolverGet(SolverSolid,LinearSolverSolid,Err)
790!   CALL CMISSSolver_LinearTypeSet(LinearSolverSolid,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
791
792
793  !Darcy
794  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
795    & SolverDarcyIndex,DynamicSolverDarcy,Err)
796  CALL CMISSSolver_OutputTypeSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE,Err)
797  CALL CMISSSolver_DynamicThetaSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_THETA,Err)
798!   CALL CMISSSolverDynamicDynamicSet(DynamicSolverDarcy,.TRUE.,Err)
799  CALL CMISSSolver_DynamicLinearSolverGet(DynamicSolverDarcy,LinearSolverDarcy,Err)
800  IF(LINEAR_SOLVER_DARCY_DIRECT_FLAG) THEN
801    CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
802    CALL CMISSSolver_LibraryTypeSet(LinearSolverDarcy,CMISS_SOLVER_MUMPS_LIBRARY,Err)
803  ELSE
804    CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
805    CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverDarcy,MAXIMUM_ITERATIONS,Err)
806    CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverDarcy,DIVERGENCE_TOLERANCE,Err)
807    CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverDarcy,RELATIVE_TOLERANCE,Err)
808    CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverDarcy,ABSOLUTE_TOLERANCE,Err)
809    CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverDarcy,RESTART_VALUE,Err)
810  ENDIF
811
812  CALL CMISSProblem_SolversCreateFinish(Problem,Err)
813
814  !
815  !================================================================================================================================
816  !
817
818  !SOLVER EQUATIONS
819
820  CALL CMISSSolver_Initialise(SolverSolid,Err)
821  CALL CMISSSolver_Initialise(LinearSolverDarcy,Err)
822
823  CALL CMISSSolverEquations_Initialise(SolverEquationsSolid,Err)
824  CALL CMISSSolverEquations_Initialise(SolverEquationsDarcy,Err)
825
826  CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
827  !
828  !Get the finite elasticity solver equations
829  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), &
830    & SolverSolidIndex,SolverSolid,Err)
831  CALL CMISSSolver_SolverEquationsGet(SolverSolid,SolverEquationsSolid,Err)
832  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsSolid,CMISS_SOLVER_SPARSE_MATRICES,Err)
833  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsSolid,EquationsSetSolid,EquationsSetIndex,Err)
834  !
835  !Get the Darcy solver equations
836  CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), &
837    & SolverDarcyIndex,LinearSolverDarcy,Err)
838  CALL CMISSSolver_SolverEquationsGet(LinearSolverDarcy,SolverEquationsDarcy,Err)
839  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsDarcy,CMISS_SOLVER_SPARSE_MATRICES,Err)
840  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsDarcy,EquationsSetDarcy,EquationsSetIndex,Err)
841  !
842  CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
843
844  !
845  !================================================================================================================================
846  !
847
848  !------------------------------------
849  ! ASSIGN BOUNDARY CONDITIONS - SOLID (absolute nodal parameters)
850  !Solid is computed in absolute position, rather than displacement. Thus BCs for absolute position
851  CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsSolid,Err)
852  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsSolid,BoundaryConditionsSolid,Err)
853
854!   !Get surfaces
855!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_FRONT_SURFACE, &
856!     & Face1Nodes,FaceXi(1),Err)
857!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BACK_SURFACE, &
858!     & Face2Nodes,FaceXi(2),Err)
859!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_RIGHT_SURFACE, &
860!     & Face3Nodes,FaceXi(3),Err)
861!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_LEFT_SURFACE, &
862!     & Face4Nodes,FaceXi(4),Err)
863!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_TOP_SURFACE, &
864!     & Face5Nodes,FaceXi(5),Err)
865!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,SolidDisplMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BOTTOM_SURFACE, &
866!     & Face6Nodes,FaceXi(6),Err)
867
868! ! write(*,'("Face1Nodes=[",100I3,"]")') Face1Nodes
869! ! write(*,'("Face2Nodes=[",100I3,"]")') Face2Nodes
870! ! write(*,'("Face3Nodes=[",100I3,"]")') Face3Nodes
871! ! write(*,'("Face4Nodes=[",100I3,"]")') Face4Nodes
872! ! write(*,'("Face5Nodes=[",100I3,"]")') Face5Nodes
873! ! write(*,'("Face6Nodes=[",100I3,"]")') Face6Nodes
874! ! write(*,'("FaceXi=[",100I3,"]")') FaceXi
875!
876allocate(Face1Nodes(21),Face2Nodes(21),Face3Nodes(21),Face4Nodes(21),Face5Nodes(9),Face6Nodes(9))
877Face1Nodes=[  1, 17,  2, 22, 23, 24,  5, 31,  6, 36, 37, 38,  9, 45, 10, 50, 51, 52, 13, 59, 14]
878Face2Nodes=[  3, 21,  4, 28, 29, 30,  7, 35,  8, 42, 43, 44, 11, 49, 12, 56, 57, 58, 15, 63, 16]
879Face3Nodes=[  2, 20,  4, 24, 27, 30,  6, 34,  8, 38, 41, 44, 10, 48, 12, 52, 55, 58, 14, 62, 16]
880Face4Nodes=[  1, 18,  3, 22, 25, 28,  5, 32,  7, 36, 39, 42,  9, 46, 11, 50, 53, 56, 13, 60, 15]
881Face5Nodes=[ 13, 59, 14, 60, 61, 62, 15, 63, 16]
882Face6Nodes=[  1, 17,  2, 18, 19, 20,  3, 21,  4]
883FaceXi=[ -2,  2,  1, -1,  3, -3]
884
885
886  ! Fix the bottom in z direction
887  DO NN=1,SIZE(Face6Nodes,1)
888    NODE=Face6Nodes(NN)
889    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
890    IF(NodeDomain==ComputationalNodeNumber) THEN
891      CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, &
892        & ZCoord,Err)
893      CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, &
894        & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err)
895      WRITE(*,*) "FIXING NODE",NODE,"AT BOTTOM IN Z DIRECTION"
896    ENDIF
897  ENDDO
898
899  ! Fix the top in z direction
900  DO NN=1,SIZE(Face5Nodes,1)
901    NODE=Face5Nodes(NN)
902    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
903    IF(NodeDomain==ComputationalNodeNumber) THEN
904      CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, &
905        & ZCoord,Err)
906      CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, &
907        & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err)
908      WRITE(*,*) "FIXING NODE",NODE,"AT TOP IN Z DIRECTION"
909    ENDIF
910  ENDDO
911
912  !Fix more nodes at the bottom to stop free body motion
913  X_FIXED=.FALSE.
914  Y_FIXED=.FALSE.
915  DO NN=1,SIZE(Face6Nodes,1)
916    NODE=Face6Nodes(NN)
917    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
918    IF(NodeDomain==ComputationalNodeNumber) THEN
919      CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,1, &
920        & XCoord,Err)
921      CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,2, &
922        & YCoord,Err)
923
924      !Fix Origin displacement in x and y (z already fixed)
925      IF(ABS(XCoord)<1.0E-6_CMISSDP.AND.ABS(YCoord)<1.0E-6_CMISSDP) THEN
926        CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,1, &
927          & CMISS_BOUNDARY_CONDITION_FIXED,XCoord,Err)
928        CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,2, &
929          & CMISS_BOUNDARY_CONDITION_FIXED,YCoord,Err)
930        WRITE(*,*) "FIXING ORIGIN NODE",NODE,"IN X AND Y DIRECTION"
931        X_FIXED=.TRUE.
932        Y_FIXED=.TRUE.
933      ENDIF
934
935      !Fix nodal displacements at (X_DIM,0) in y
936      IF(ABS(XCoord - X_DIM)<1.0E-6_CMISSDP .AND. ABS(YCoord)<1.0E-6_CMISSDP) THEN
937        CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,2, &
938          & CMISS_BOUNDARY_CONDITION_FIXED,YCoord,Err)
939        WRITE(*,*) "FIXING NODES",NODE,"AT (X_DIM,0) IN Y DIRECTION"
940        Y_FIXED=.TRUE.
941      ENDIF
942
943      !Fix nodal displacements at (0,Y_DIM) in x
944      IF(ABS(XCoord)<1.0E-6_CMISSDP .AND. ABS(YCoord - Y_DIM)<1.0E-6_CMISSDP) THEN
945        CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsSolid,DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,1, &
946          & CMISS_BOUNDARY_CONDITION_FIXED,XCoord,Err)
947        WRITE(*,*) "FIXING NODES",NODE,"AT (0,Y_DIM) IN X DIRECTION"
948        X_FIXED=.TRUE.
949      ENDIF
950
951    ENDIF
952  ENDDO
953!   CALL MPI_REDUCE(X_FIXED,X_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
954!   CALL MPI_REDUCE(Y_FIXED,Y_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
955!   IF(ComputationalNodeNumber==0) THEN
956!     IF(.NOT.(X_OKAY.AND.Y_OKAY)) THEN
957!       WRITE(*,*) "Free body motion could not be prevented!"
958!       CALL CMISSFinalise(Err)
959!       STOP
960!     ENDIF
961!   ENDIF
962
963  CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsSolid,Err)
964  !------------------------------------
965
966
967  !------------------------------------
968  ! ASSIGN BOUNDARY CONDITIONS - FLUID
969  CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsDarcy,Err)
970  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsDarcy,BoundaryConditionsDarcy,Err)
971
972!   !Get surfaces
973!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_FRONT_SURFACE, &
974!     & Face7Nodes,FaceXi(1),Err)
975!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BACK_SURFACE, &
976!     & Face8Nodes,FaceXi(2),Err)
977!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_RIGHT_SURFACE, &
978!     & Face9Nodes,FaceXi(3),Err)
979!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_LEFT_SURFACE, &
980!     & Face10Nodes,FaceXi(4),Err)
981!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_TOP_SURFACE, &
982!     & Face11Nodes,FaceXi(5),Err)
983!   CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,DarcyVelMeshComponentNumber,CMISS_GENERATED_MESH_REGULAR_BOTTOM_SURFACE, &
984!     & Face12Nodes,FaceXi(6),Err)
985
986! ! write(*,'("Face7Nodes=[",100I3,"]")') Face7Nodes
987! ! write(*,'("Face8Nodes=[",100I3,"]")') Face8Nodes
988! ! write(*,'("Face9Nodes=[",100I3,"]")') Face9Nodes
989!

Large files files are truncated, but you can click here to view the full file