/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenMultiCompDarcy/src/QuadraticEllipsoidDrivenMultiCompDarcyExample.f90
FORTRAN Modern | 1318 lines | 824 code | 154 blank | 340 comment | 0 complexity | 61ea8a05237d93b94245a785e439c9ef MD5 | raw file
Large files files are truncated, but you can click here to view the full file
1!> \file 2!> \author Christian Michler 3!> \brief This is an example program to solve a coupled Finite Elastiticity Darcy equation on a cylindrical geometry 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): Christian Michler, Jack Lee 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/QuadraticEllipsoidDrivenDarcy/src/QuadraticEllipsoidDrivenDarcyExample.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/QuadraticEllipsoidDrivenDarcy/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/QuadraticEllipsoidDrivenDarcy/build-intel'>Linux GNU Build</a> 47!! 48!< 49 50! ! 51! ! This example considers a coupled Finite Elasticity Darcy problem on a cylindrical geometry 52! ! 53 54!> Main program 55 56PROGRAM QUADRATICELLIPSOIDDRIVENMULTICOMPDARCYEXAMPLE 57 58 ! 59 !================================================================================================================================ 60 ! 61 62 !PROGRAM LIBRARIES 63 64 USE OPENCMISS 65 USE MPI 66 67#ifdef WIN32 68 USE IFQWINCMISS 69#endif 70 71 ! 72 !================================================================================================================================ 73 ! 74 75 !PROGRAM VARIABLES AND TYPES 76 77 IMPLICIT NONE 78 79 !Test program parameters 80 81 REAL(CMISSDP), PARAMETER :: PI=3.14159_CMISSDP 82 REAL(CMISSDP), PARAMETER :: LONG_AXIS=2.0_CMISSDP 83 REAL(CMISSDP), PARAMETER :: SHORT_AXIS=1.0_CMISSDP 84 REAL(CMISSDP), PARAMETER :: WALL_THICKNESS=0.5_CMISSDP 85 REAL(CMISSDP), PARAMETER :: CUTOFF_ANGLE=1.5708_CMISSDP 86 REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_INTERSECTION=1.73205_CMISSDP !Slope of fibres in base endocardium = 60 degrees 87 REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_CHANGE=-3.4641_CMISSDP !Slope change of fibres from 60 to -60 degrees in transmural direction 88 REAL(CMISSDP), PARAMETER :: SHEET_SLOPE_BASE_ENDO=1.0_CMISSDP !Slope of sheet at base endocardium 89 90 REAL(CMISSDP), PARAMETER :: INNER_PRESSURE=2.0_CMISSDP !Positive is compressive 91 REAL(CMISSDP), PARAMETER :: OUTER_PRESSURE=0.0_CMISSDP !Positive is compressive 92 REAL(CMISSDP), PARAMETER :: C1= 2.0_CMISSDP 93 REAL(CMISSDP), PARAMETER :: C2= 6.0_CMISSDP 94 REAL(CMISSDP), PARAMETER :: C3=10.0_CMISSDP 95 96 INTEGER(CMISSIntg), PARAMETER :: NumberGlobalXElements=4 ! X ==NUMBER_GLOBAL_CIRCUMFERENTIAL_ELEMENTS 97 INTEGER(CMISSIntg), PARAMETER :: NumberGlobalYElements=4 ! Y ==NUMBER_GLOBAL_LONGITUDINAL_ELEMENTS 98 INTEGER(CMISSIntg), PARAMETER :: NumberGlobalZElements=1 ! Z ==NUMBER_GLOBAL_TRANSMURAL_ELEMENTS 99 INTEGER(CMISSIntg) :: NumberOfDomains 100 101 INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1 102 INTEGER(CMISSIntg), PARAMETER :: NumberOfSpatialCoordinates=3 103 INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=1 104 INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=1 105 INTEGER(CMISSIntg), PARAMETER :: QuadraticCollapsedBasisUserNumber=2 106 INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=3 107 INTEGER(CMISSIntg), PARAMETER :: LinearCollapsedBasisUserNumber=4 108 INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=1 109 INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=2 110 INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=1 111 INTEGER(CMISSIntg), PARAMETER :: DerivativeUserNumber=1 112 113 INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshDimensions=3 114 INTEGER(CMISSIntg), PARAMETER :: NumberOfXiCoordinates=3 115 INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshComponents=2 116 INTEGER(CMISSIntg), PARAMETER :: QuadraticMeshComponentNumber=1 117 INTEGER(CMISSIntg), PARAMETER :: LinearMeshComponentNumber=2 118 INTEGER(CMISSIntg), PARAMETER :: TotalNumberOfElements=1 119 120 INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumberSolid=1 121 INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumberDarcy=11 122 INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1 123 INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3 124 125 INTEGER(CMISSIntg), PARAMETER :: FieldFibreUserNumber=2 126 INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfVariables=1 127 INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfComponents=3 128 129 INTEGER(CMISSIntg), PARAMETER :: FieldMaterialUserNumber=3 130 INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfVariables=1 131 INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfComponents=3 !2 132 133 INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidUserNumber=4 134 INTEGER(CMISSIntg) :: FieldDependentSolidNumberOfVariables 135 INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4 136 INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4 137 138 INTEGER(CMISSIntg), PARAMETER :: EquationSetSolidUserNumber=55 139 INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldSolidUserNumber=25 140 INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=1 141 INTEGER(CMISSIntg) :: EquationsSetFieldUserNumberDarcy 142 INTEGER(CMISSIntg) :: icompartment,Ncompartments,num_var,componentnum,Nparams 143 144 INTEGER(CMISSIntg) :: MaterialsFieldUserNumberDarcy 145 INTEGER(CMISSIntg) :: EquationsSetUserNumberDarcy 146 147 148 INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1 149 INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2 150 INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1 151 INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1 152 INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1 153 INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPorosity=1 154 INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcyPermOverVis=2 155 INTEGER(CMISSIntg), PARAMETER :: SolidDisplMeshComponentNumber=1 156 INTEGER(CMISSIntg), PARAMETER :: SolidLagrMultMeshComponentNumber=2 157 INTEGER(CMISSIntg), PARAMETER :: SolidGeometryMeshComponentNumber=SolidDisplMeshComponentNumber 158 !Program types 159 160 161 !Program variables 162 163 INTEGER(CMISSIntg) :: MPI_IERROR 164 INTEGER(CMISSIntg) :: EquationsSetIndex 165 INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber 166 REAL(CMISSDP) :: FibreFieldAngle(3) 167 REAL(CMISSDP) :: nu,theta,omega,XI3,XI3delta,XI2delta, zero 168 INTEGER(CMISSIntg) ::i,j,k,component_idx,node_idx,TOTAL_NUMBER_NODES_XI(3) 169 !For grabbing surfaces 170 INTEGER(CMISSIntg) :: InnerNormalXi,OuterNormalXi,TopNormalXi 171 INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodes(:) 172 INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodes(:) 173 INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodes(:) 174 INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodesDarcyVel(:) 175 INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodesDarcyVel(:) 176 INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodesDarcyVel(:) 177 INTEGER(CMISSIntg) :: NN,NODE,NodeDomain 178 REAL(CMISSDP) :: XCoord,YCoord,ZCoord 179 LOGICAL :: X_FIXED,Y_FIXED,X_OKAY,Y_OKAY 180 181 INTEGER(CMISSIntg) :: GeometricFieldDarcyMeshComponentNumber, DarcyVelMeshComponentNumber, DarcyMassIncreaseMeshComponentNumber 182 INTEGER(CMISSIntg) :: MeshComponentNumber_dummy 183 184 INTEGER(CMISSIntg) :: NUMBER_OF_DOMAINS 185 186 INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS 187 188 INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS 189 INTEGER(CMISSIntg) :: RESTART_VALUE 190 191 INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT 192 INTEGER(CMISSIntg) :: COMPONENT_NUMBER 193 INTEGER(CMISSIntg) :: CONDITION 194 195 INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY 196 INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE 197 INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE 198 199 REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z 200 REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2 201 REAL(CMISSDP) :: GEOMETRY_TOLERANCE 202 INTEGER(CMISSIntg) :: EDGE_COUNT 203 INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID 204 REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4) 205 REAL(CMISSDP) :: DIVERGENCE_TOLERANCE 206 REAL(CMISSDP) :: RELATIVE_TOLERANCE 207 REAL(CMISSDP) :: ABSOLUTE_TOLERANCE 208 REAL(CMISSDP) :: LINESEARCH_ALPHA 209 REAL(CMISSDP) :: VALUE 210 REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY 211 212 LOGICAL :: EXPORT_FIELD_IO 213 LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG 214 215 !CMISS variables 216 217 TYPE(CMISSBasisType) :: QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis 218 TYPE(CMISSBoundaryConditionsType) :: BoundaryConditions 219 TYPE(CMISSCoordinateSystemType) :: CoordinateSystem, WorldCoordinateSystem 220 TYPE(CMISSMeshType) :: Mesh 221 TYPE(CMISSGeneratedMeshType) :: GeneratedMesh 222 TYPE(CMISSDecompositionType) :: Decomposition 223 TYPE(CMISSEquationsType) :: Equations 224 TYPE(CMISSEquationsSetType) :: EquationsSetSolid 225 TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid 226 TYPE(CMISSFieldType) :: DependentFieldSolid,EquationsSetFieldSolid 227 228 TYPE(CMISSFieldType) :: GeometricFieldDarcy 229 230 TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: MaterialsFieldDarcy 231 TYPE(CMISSFieldType), ALLOCATABLE, DIMENSION(:) :: EquationsSetFieldDarcy 232 !Boundary conditions 233 TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy 234 !Equations sets 235 TYPE(CMISSEquationsSetType), ALLOCATABLE, DIMENSION(:) :: EquationsSetDarcy 236 !Equations 237 TYPE(CMISSEquationsType), ALLOCATABLE, DIMENSION(:) :: EquationsDarcy 238 239 240 241 TYPE(CMISSFieldsType) :: Fields 242 TYPE(CMISSProblemType) :: Problem 243 TYPE(CMISSRegionType) :: Region,WorldRegion 244 TYPE(CMISSSolverType) :: SolverSolid,LinearSolverSolid 245 TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid 246 TYPE(CMISSControlLoopType) :: ControlLoop 247 248 !Solvers 249 TYPE(CMISSSolverType) :: DynamicSolverDarcy 250 TYPE(CMISSSolverType) :: LinearSolverDarcy 251 !Solver equations 252 TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy 253 254 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME 255 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME 256 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA 257 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT 258 259#ifdef WIN32 260 !Quickwin type 261 LOGICAL :: QUICKWIN_STATUS=.FALSE. 262 TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG 263#endif 264 265 !Generic CMISS variables 266 INTEGER(CMISSIntg) :: Err 267 INTEGER(CMISSIntg), ALLOCATABLE, DIMENSION(:) :: VariableTypes 268 REAL(CMISSDP), ALLOCATABLE, DIMENSION(:,:) :: CouplingCoeffs,ConstitutiveParams 269#ifdef WIN32 270 !Initialise QuickWin 271 QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title 272 QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows 273 QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN 274 !Set the window parameters 275 QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG) 276 !If attempt fails set with system estimated values 277 IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG) 278#endif 279 280 !Set initial values 281 NUMBER_OF_DIMENSIONS=3_CMISSIntg 282 INITIAL_FIELD_DARCY(1)=0.0_CMISSDP 283 INITIAL_FIELD_DARCY(2)=0.0_CMISSDP 284 INITIAL_FIELD_DARCY(3)=0.0_CMISSDP 285 INITIAL_FIELD_DARCY(4)=0.0_CMISSDP 286 !Set material parameters 287 POROSITY_PARAM_DARCY=0.1_CMISSDP 288 PERM_OVER_VIS_PARAM_DARCY=1.0_CMISSDP 289 !Set output parameter 290 !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput) 291 DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_PROGRESS_OUTPUT 292 LINEAR_SOLVER_DARCY_OUTPUT_TYPE=CMISS_SOLVER_SOLVER_OUTPUT 293 !(NoOutput/TimingOutput/MatrixOutput/ElementOutput) 294 EQUATIONS_DARCY_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT 295 296 !Set time parameter 297 DYNAMIC_SOLVER_DARCY_START_TIME=0.0_CMISSDP 298 DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP 299 DYNAMIC_SOLVER_DARCY_STOP_TIME=5_CMISSIntg * DYNAMIC_SOLVER_DARCY_TIME_INCREMENT 300 DYNAMIC_SOLVER_DARCY_THETA=1.0_CMISSDP !2.0_CMISSDP/3.0_CMISSDP 301 !Set result output parameter 302 DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY=1 303 !Set solver parameters 304 LINEAR_SOLVER_DARCY_DIRECT_FLAG=.TRUE. 305 RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP 306 ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP 307 DIVERGENCE_TOLERANCE=1.0E5_CMISSDP !default: 1.0E5 308 MAXIMUM_ITERATIONS=10000_CMISSIntg !default: 100000 309 RESTART_VALUE=30_CMISSIntg !default: 30 310 LINESEARCH_ALPHA=1.0_CMISSDP 311 icompartment =1_CMISSIntg 312 Ncompartments=4_CMISSIntg 313 314 !LinearMeshComponentNumber/QuadraticMeshComponentNumber 315 DarcyVelMeshComponentNumber = LinearMeshComponentNumber 316 DarcyMassIncreaseMeshComponentNumber = LinearMeshComponentNumber 317! GeometricFieldDarcyMeshComponentNumber = DarcyVelMeshComponentNumber 318 GeometricFieldDarcyMeshComponentNumber = QuadraticMeshComponentNumber 319 320 ALLOCATE (EquationsSetDarcy(Ncompartments)) 321 ALLOCATE (EquationsSetFieldDarcy(Ncompartments)) 322 ALLOCATE (MaterialsFieldDarcy(Ncompartments)) 323 ALLOCATE (EquationsDarcy(Ncompartments)) 324 ! 325 !================================================================================================================================ 326 ! 327 328 !Intialise cmiss 329 CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err) 330 331 CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err) 332 333 WRITE(*,'(A)') "Program starting." 334 335 !Set all diganostic levels on for testing 336 CALL CMISSDiagnosticsSetOn(CMISS_FROM_DIAG_TYPE,(/1,2,3,4,5/),"Diagnostics",(/"PROBLEM_FINITE_ELEMENT_CALCULATE"/),Err) 337 338 !Get the number of computational nodes and this computational node number 339 CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err) 340 CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err) 341 342 !Broadcast the number of elements in the X,Y and Z directions and the number of partitions to the other computational nodes 343! CALL MPI_BCAST(NumberGlobalXElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 344! CALL MPI_BCAST(NumberGlobalYElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 345! CALL MPI_BCAST(NumberGlobalZElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 346! CALL MPI_BCAST(NumberOfDomains,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 347 NumberOfDomains=NumberOfComputationalNodes 348 349 !Create a CS - default is 3D rectangular cartesian CS with 0,0,0 as origin 350 CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err) 351 CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err) 352 CALL CMISSCoordinateSystem_TypeSet(CoordinateSystem,CMISS_COORDINATE_RECTANGULAR_CARTESIAN_TYPE,Err) 353 CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NumberOfSpatialCoordinates,Err) 354 CALL CMISSCoordinateSystem_OriginSet(CoordinateSystem,(/0.0_CMISSDP,0.0_CMISSDP,0.0_CMISSDP/),Err) 355 CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err) 356 357 ! 358 !================================================================================================================================ 359 ! 360 361 !Create a region and assign the CS to the region 362 CALL CMISSRegion_Initialise(Region,Err) 363 CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err) 364 CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err) 365 CALL CMISSRegion_CreateFinish(Region,Err) 366 367 ! 368 !================================================================================================================================ 369 ! 370 371 !Define basis functions - tri-linear Lagrange and tri-Quadratic Lagrange, each with collapsed variant 372 !Quadratic Basis 373 CALL CMISSBasis_Initialise(QuadraticBasis,Err) 374 CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err) 375 CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, & 376 & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err) 377 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, & 378 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 379! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 380 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err) !Have to do this 381 CALL CMISSBasis_CreateFinish(QuadraticBasis,Err) 382 383 !Collapsed Quadratic Basis 384 CALL CMISSBasis_Initialise(QuadraticCollapsedBasis,Err) 385 CALL CMISSBasis_CreateStart(QuadraticCollapsedBasisUserNumber,QuadraticCollapsedBasis,Err) 386 CALL CMISSBasis_TypeSet(QuadraticCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err) 387 CALL CMISSBasis_NumberOfXiSet(QuadraticCollapsedBasis,NumberOfXiCoordinates,Err) 388 CALL CMISSBasis_InterpolationXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, & 389 & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err) 390 CALL CMISSBasis_CollapsedXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED, & 391 & CMISS_BASIS_COLLAPSED_AT_XI0,CMISS_BASIS_NOT_COLLAPSED/),Err) 392 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticCollapsedBasis, & 393 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 394! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 395 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticCollapsedBasis,.true.,Err) !Have to do this 396 CALL CMISSBasis_CreateFinish(QuadraticCollapsedBasis,Err) 397 398 !Linear Basis 399 CALL CMISSBasis_Initialise(LinearBasis,Err) 400 CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err) 401 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, & 402 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 403! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 404 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err) !Have to do this (unused) due to field_interp setup 405 CALL CMISSBasis_CreateFinish(LinearBasis,Err) 406 407 !Collapsed Linear Basis 408 CALL CMISSBasis_Initialise(LinearCollapsedBasis,Err) 409 CALL CMISSBasis_CreateStart(LinearCollapsedBasisUserNumber,LinearCollapsedBasis,Err) 410 CALL CMISSBasis_TypeSet(LinearCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err) 411 CALL CMISSBasis_NumberOfXiSet(LinearCollapsedBasis,NumberOfXiCoordinates,Err) 412 CALL CMISSBasis_InterpolationXiSet(LinearCollapsedBasis,(/CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION, & 413 & CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION,CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION/),Err) 414 CALL CMISSBasis_CollapsedXiSet(LinearCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED,CMISS_BASIS_COLLAPSED_AT_XI0, & 415 & CMISS_BASIS_NOT_COLLAPSED/),Err) 416 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearCollapsedBasis, & 417 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 418! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 419 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearCollapsedBasis,.true.,Err) !Have to do this (unused) due to field_interp setup 420 CALL CMISSBasis_CreateFinish(LinearCollapsedBasis,Err) 421 422 ! 423 !================================================================================================================================ 424 ! 425 426 !Start the creation of a generated ellipsoid mesh 427 CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err) 428 CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err) 429 !Set up an ellipsoid mesh 430 CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_MESH_TYPE,Err) 431 !Set the quadratic and linear bases 432 CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,[QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis],Err) 433 !Define the mesh on the region 434 CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/LONG_AXIS,SHORT_AXIS,WALL_THICKNESS,CUTOFF_ANGLE/),Err) 435 CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NumberGlobalXElements,NumberGlobalYElements, & 436 & NumberGlobalZElements/),Err) 437 438 !Finish the creation of a generated mesh in the region 439 CALL CMISSMesh_Initialise(Mesh,Err) 440 CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err) 441 442 ! 443 !================================================================================================================================ 444 ! 445 446 !Create a decomposition 447 CALL CMISSDecomposition_Initialise(Decomposition,Err) 448 CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err) 449 CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err) 450 CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err) 451 CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.TRUE.,Err) 452 CALL CMISSDecomposition_CreateFinish(Decomposition,Err) 453 454 ! 455 !================================================================================================================================ 456 ! 457 458 ! --- GeometricFieldSolid --- 459 !Create a field to put the geometry (default is geometry) 460 CALL CMISSField_Initialise(GeometricFieldSolid,Err) 461 CALL CMISSField_CreateStart(FieldGeometryUserNumberSolid,Region,GeometricFieldSolid,Err) 462 CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err) 463 CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err) 464 CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometryNumberOfVariables,Err) 465 CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err) 466 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err) 467 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err) 468 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err) 469 CALL CMISSField_CreateFinish(GeometricFieldSolid,Err) 470 471 !Update the geometric field parameters 472 CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err) 473 474 ! 475 !================================================================================================================================ 476 ! 477 478 ! --- GeometricFieldDarcy --- 479 !Create a field to put the geometry (default is geometry) 480 CALL CMISSField_Initialise(GeometricFieldDarcy,Err) 481 CALL CMISSField_CreateStart(FieldGeometryUserNumberDarcy,Region,GeometricFieldDarcy,Err) 482 CALL CMISSField_MeshDecompositionSet(GeometricFieldDarcy,Decomposition,Err) 483 CALL CMISSField_TypeSet(GeometricFieldDarcy,CMISS_FIELD_GEOMETRIC_TYPE,Err) 484 CALL CMISSField_NumberOfVariablesSet(GeometricFieldDarcy,FieldGeometryNumberOfVariables,Err) 485 CALL CMISSField_NumberOfComponentsSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err) 486 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,1, & 487 & GeometricFieldDarcyMeshComponentNumber,Err) 488 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,2, & 489 & GeometricFieldDarcyMeshComponentNumber,Err) 490 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,3, & 491 & GeometricFieldDarcyMeshComponentNumber,Err) 492 CALL CMISSField_CreateFinish(GeometricFieldDarcy,Err) 493 494 !Update the geometric field parameters 495 CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldDarcy,Err) 496 497 ! 498 !================================================================================================================================ 499 ! 500 501 !Create a fibre field and attach it to the geometric field 502 CALL CMISSField_Initialise(FibreFieldSolid,Err) 503 CALL CMISSField_CreateStart(FieldFibreUserNumber,Region,FibreFieldSolid,Err) 504 CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err) 505 CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err) 506 CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err) 507 CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreNumberOfVariables,Err) 508 CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreNumberOfComponents,Err) 509 CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err) 510 CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err) 511 CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err) 512 CALL CMISSField_CreateFinish(FibreFieldSolid,Err) 513 514 !Set Fibre directions (this block is parallel-untested) 515 node_idx=0 516 !This is valid only for quadratic basis functions 517 TOTAL_NUMBER_NODES_XI(1)=NumberGlobalXElements*2 518 TOTAL_NUMBER_NODES_XI(2)=NumberGlobalYElements*2+1 519 TOTAL_NUMBER_NODES_XI(3)=NumberGlobalZElements*2+1 520 521 XI2delta=(PI-CUTOFF_ANGLE)/(TOTAL_NUMBER_NODES_XI(2)-1) 522 XI3=0 523 XI3delta=(1.0)/(TOTAL_NUMBER_NODES_XI(3)-1) 524 zero=0 525 DO k=1, TOTAL_NUMBER_NODES_XI(3) 526 !Apex nodes 527 j=1 528 i=1 529 node_idx=node_idx+1 530 CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err) 531 IF(NodeDomain==ComputationalNodeNumber) THEN 532 FibreFieldAngle=(/zero,zero,zero/) 533 DO component_idx=1,FieldFibreNumberOfComponents 534 CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, & 535 & DerivativeUserNumber,node_idx,component_idx,FibreFieldAngle(component_idx),Err) 536 ENDDO 537 ENDIF 538 theta=atan(FIBRE_SLOPE_CHANGE*XI3+FIBRE_SLOPE_INTERSECTION) 539 DO j=2, TOTAL_NUMBER_NODES_XI(2) 540 nu=PI-XI2delta*(j-1) 541 omega=PI/2+cos(2*nu)*atan(SHEET_SLOPE_BASE_ENDO)*(-2*XI3+1) 542 DO i=1, TOTAL_NUMBER_NODES_XI(1) 543 node_idx=node_idx+1 544 CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err) 545 IF(NodeDomain==ComputationalNodeNumber) THEN 546 FibreFieldAngle=(/theta,zero,omega/) 547 DO component_idx=1,FieldFibreNumberOfComponents 548 CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, & 549 & DerivativeUserNumber, node_idx,component_idx,FibreFieldAngle(component_idx),Err) 550 ENDDO 551 ENDIF 552 ENDDO 553 ENDDO 554 XI3=XI3+XI3delta 555 ENDDO 556 557 !Create a material field and attach it to the geometric field 558 CALL CMISSField_Initialise(MaterialFieldSolid,Err) 559 CALL CMISSField_CreateStart(FieldMaterialUserNumber,Region,MaterialFieldSolid,Err) 560 CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err) 561 CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err) 562 CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err) 563 CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialNumberOfVariables,Err) 564 CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialNumberOfComponents,Err) 565 CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_CONSTANT_INTERPOLATION,Err) 566 CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_CONSTANT_INTERPOLATION,Err) 567 CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_CONSTANT_INTERPOLATION,Err) 568 CALL CMISSField_CreateFinish(MaterialFieldSolid,Err) 569 570 !Set Mooney-Rivlin constants 571 CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,C1,Err) 572 CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,C2,Err) 573 CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,C3,Err) 574 575 ! 576 !================================================================================================================================ 577 ! 578 579 !EQUATIONS SETS 580 DO icompartment = 1,Ncompartments 581 EquationsSetFieldUserNumberDarcy = 100_CMISSIntg+icompartment 582 EquationsSetUserNumberDarcy = 200_CMISSIntg+icompartment 583 !Create the equations set for ALE Darcy 584 CALL CMISSField_Initialise(EquationsSetFieldDarcy(icompartment),Err) 585 CALL CMISSEquationsSet_Initialise(EquationsSetDarcy(icompartment),Err) 586 CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricFieldDarcy, & 587 & CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, & 588 & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,& 589 & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy(icompartment),EquationsSetDarcy(icompartment),Err) 590 !Finish creating the equations set 591 CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy(icompartment),Err) 592 !Set the values for the equations set field to be the current compartment number (1 - N), and the total number of compartments (N) 593 CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, & 594 & CMISS_FIELD_VALUES_SET_TYPE,1,icompartment,Err) 595 CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, & 596 & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err) 597 ENDDO 598 599 !-------------------------------------------------------------------------------------------------------------------------------- 600 ! Solid 601 602 !Create the equations_set 603 CALL CMISSField_Initialise(EquationsSetFieldSolid,Err) 604 CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err) 605 CALL CMISSEquationsSet_CreateStart(EquationSetSolidUserNumber,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, & 606 & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,& 607 & EquationsSetFieldSolidUserNumber,EquationsSetFieldSolid,EquationsSetSolid,Err) 608 CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err) 609 !Set the values for the equations set field to be the current compartment number (O for the finite elasticity equations_set), and the total number of compartments (N) 610 !Need to store number of compartments, as finite elasticity uses this to calculate the total mass increase for the constiutive law 611 CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, & 612 & CMISS_FIELD_VALUES_SET_TYPE,1,0_CMISSIntg,Err) 613 CALL CMISSField_ParameterSetUpdateConstant(EquationsSetFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, & 614 & CMISS_FIELD_VALUES_SET_TYPE,2,Ncompartments,Err) 615 616 ! 617 !================================================================================================================================ 618 ! 619 ! Solid 620 621 !Create a dependent field with two variables and four components 622 CALL CMISSField_Initialise(DependentFieldSolid,Err) 623 ! 624 CALL CMISSField_CreateStart(FieldDependentSolidUserNumber,Region,DependentFieldSolid,Err) 625 ! 626 CALL CMISSField_TypeSet(DependentFieldSolid,CMISS_FIELD_GENERAL_TYPE,Err) 627 CALL CMISSField_MeshDecompositionSet(DependentFieldSolid,Decomposition,Err) 628 CALL CMISSField_GeometricFieldSet(DependentFieldSolid,GeometricFieldSolid,Err) 629 CALL CMISSField_DependentTypeSet(DependentFieldSolid,CMISS_FIELD_DEPENDENT_TYPE,Err) 630 !Create 2N+2 number of variables - 2 for solid, 2N for N Darcy compartments 631 FieldDependentSolidNumberOfVariables=2*Ncompartments+2 632 CALL CMISSField_NumberOfVariablesSet(DependentFieldSolid,FieldDependentSolidNumberOfVariables,Err) 633 !create two variables for each compartment 634 ALLOCATE(VariableTypes(2*Ncompartments+2)) 635 DO num_var=1,Ncompartments+1 636 VariableTypes(2*num_var-1)=CMISS_FIELD_U_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1)) 637 VariableTypes(2*num_var)=CMISS_FIELD_DELUDELN_VARIABLE_TYPE+(CMISS_FIELD_NUMBER_OF_VARIABLE_SUBTYPES*(num_var-1)) 638 ENDDO 639 CALL CMISSField_VariableTypesSet(DependentFieldSolid,VariableTypes,Err) 640! CALL CMISSField_VariableTypesSet(DependentFieldSolid,(/CMISS_FIELD_U_VARIABLE_TYPE, & 641! & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err) 642 CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, & 643 & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err) 644 CALL CMISSField_DimensionSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, & 645 & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err) 646 CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err) 647 CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, & 648 & FieldDependentSolidNumberOfComponents,Err) 649 DO icompartment=3,2*Ncompartments+2 650 CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err) 651 ENDDO 652! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 653! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 654! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 655 656 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,SolidDisplMeshComponentNumber,Err) 657 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,SolidDisplMeshComponentNumber,Err) 658 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,SolidDisplMeshComponentNumber,Err) 659 CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4, & 660 & CMISS_FIELD_NODE_BASED_INTERPOLATION, & 661 & Err) 662! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err) 663! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err) 664 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,4,SolidLagrMultMeshComponentNumber,Err) 665! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, & 666! & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 667! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, & 668! & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 669! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, & 670! & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 671 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1, & 672 & SolidDisplMeshComponentNumber, & 673 & Err) 674 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2, & 675 & SolidDisplMeshComponentNumber, & 676 & Err) 677 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3, & 678 & SolidDisplMeshComponentNumber, & 679 & Err) 680 CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, & 681 & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 682! CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, & 683! & CMISS_FIELD_ELEMENT_BASED_INTERPOLATION,Err) 684! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,SolidMeshComponenetNumber,Err) 685 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4, & 686 & SolidLagrMultMeshComponentNumber, & 687 & Err) 688 !loop over the number of compartments 689 DO icompartment=3,2*Ncompartments+2 690! CALL CMISSField_DimensionSet(DependentFieldSolid,VariableTypes(icompartment), & 691! & CMISS_FIELD_VECTOR_DIMENSION_TYPE,Err) 692 !CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,VariableTypes(icompartment),FieldDependentFluidNumberOfComponents,Err) 693 DO componentnum=1,FieldDependentFluidNumberOfComponents-1 694 !set dimension type 695! CALL CMISSField_DimensionSet(DependentField,VariableTypes(icompartment), & 696! & CMISS_FIELD_SCALAR_DIMENSION_TYPE,Err) 697 CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, & 698 & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 699 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment),componentnum, & 700 & DarcyVelMeshComponentNumber,Err) 701 ENDDO 702 CALL CMISSField_ComponentInterpolationSet(DependentFieldSolid,VariableTypes(icompartment), & 703 & FieldDependentFluidNumberOfComponents, & 704 & CMISS_FIELD_NODE_BASED_INTERPOLATION,Err) 705! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), & 706! & FieldDependentFluidNumberOfComponents,MESH_COMPONENT_NUMBER_PRESSURE,Err) 707 CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,VariableTypes(icompartment), & 708 & FieldDependentFluidNumberOfComponents,DarcyMassIncreaseMeshComponentNumber,Err) 709 710 ENDDO 711 712! CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err) 713! CALL CMISSField_NumberOfComponentsSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err) 714! !For this equation type, MESH_COMPONENT_NUMBER_PRESSURE is actually the mass increase component as the pressure is taken from the solid equations 715! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err) 716! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err) 717! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err) 718! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err) 719! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,MESH_COMPONENT_NUMBER_VELOCITY,Err) 720! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,MESH_COMPONENT_NUMBER_VELOCITY,Err) 721! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,MESH_COMPONENT_NUMBER_VELOCITY,Err) 722! CALL CMISSField_ComponentMeshComponentSet(DependentFieldSolid,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4,MESH_COMPONENT_NUMBER_PRESSURE,Err) 723 724 ! 725 CALL CMISSField_ScalingTypeSet(DependentFieldSolid,CMISS_FIELD_UNIT_SCALING,Err) 726 727 CALL CMISSField_CreateFinish(DependentFieldSolid,Err) 728 729 ! 730 !================================================================================================================================ 731 ! 732 733 CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentSolidUserNumber,DependentFieldSolid,Err) 734 CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err) 735 736 CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialUserNumber,MaterialFieldSolid,Err) 737 CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err) 738 739 ! 740 !================================================================================================================================ 741 ! 742 DO icompartment = 1,Ncompartments 743 CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy(icompartment),FieldDependentSolidUserNumber,& 744 & DependentFieldSolid,Err) 745 CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy(icompartment),Err) 746 ENDDO 747 748 DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS+1 749 CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 750 & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err) 751 CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U1_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 752 & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err) 753 CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U2_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 754 & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err) 755 CALL CMISSField_ComponentValuesInitialise(DependentFieldSolid,CMISS_FIELD_U3_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 756 & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err) 757 ENDDO 758 ! 759 !================================================================================================================================ 760 ! 761 ALLOCATE(CouplingCoeffs(Ncompartments,Ncompartments)) 762 IF(Ncompartments==2)THEN 763 CouplingCoeffs(1,1)=0.0E-01_CMISSDP 764! CouplingCoeffs(1,2)=-1.0E-04_CMISSDP 765! CouplingCoeffs(2,1)=-1.0E-04_CMISSDP 766 CouplingCoeffs(1,2)=0.0E-01_CMISSDP 767 CouplingCoeffs(2,1)=0.0E-01_CMISSDP 768 CouplingCoeffs(2,2)=0.0E-01_CMISSDP 769 ELSE IF(Ncompartments==3)THEN 770 CouplingCoeffs(1,1)=1.0E-02_CMISSDP 771 CouplingCoeffs(1,2)=1.0E-02_CMISSDP 772 CouplingCoeffs(1,3)=0.0E-02_CMISSDP 773 CouplingCoeffs(2,1)=1.0E-02_CMISSDP 774 CouplingCoeffs(2,2)=2.0E-02_CMISSDP 775 CouplingCoeffs(2,3)=1.0E-02_CMISSDP 776 CouplingCoeffs(3,1)=0.0E-02_CMISSDP 777 CouplingCoeffs(3,2)=1.0E-02_CMISSDP 778 CouplingCoeffs(3,3)=1.0E-02_CMISSDP 779 ELSE IF(Ncompartments==4)THEN 780 CouplingCoeffs(1,1)=0.0E-02_CMISSDP 781 CouplingCoeffs(1,2)=0.0E-02_CMISSDP 782 CouplingCoeffs(1,3)=0.0E-02_CMISSDP 783 CouplingCoeffs(1,4)=0.0E-02_CMISSDP 784 CouplingCoeffs(2,1)=0.0E-02_CMISSDP 785 CouplingCoeffs(2,2)=0.0E-02_CMISSDP 786 CouplingCoeffs(2,3)=0.0E-02_CMISSDP 787 CouplingCoeffs(2,4)=0.0E-02_CMISSDP 788 CouplingCoeffs(3,1)=0.0E-02_CMISSDP 789 CouplingCoeffs(3,2)=0.0E-02_CMISSDP 790 CouplingCoeffs(3,3)=0.0E-02_CMISSDP 791 CouplingCoeffs(3,4)=0.0E-02_CMISSDP 792 CouplingCoeffs(4,1)=0.0E-02_CMISSDP 793 CouplingCoeffs(4,2)=0.0E-02_CMISSDP 794 CouplingCoeffs(4,3)=0.0E-02_CMISSDP 795 CouplingCoeffs(4,4)=0.0E-02_CMISSDP 796 ELSE IF(Ncompartments==5)THEN 797 CouplingCoeffs(1,1)=0.0E-02_CMISSDP 798 CouplingCoeffs(1,2)=0.0E-02_CMISSDP 799 CouplingCoeffs(1,3)=0.0E-02_CMISSDP 800 CouplingCoeffs(1,4)=0.0E-02_CMISSDP 801 CouplingCoeffs(1,5)=0.0E-02_CMISSDP 802 CouplingCoeffs(2,1)=0.0E-02_CMISSDP 803 CouplingCoeffs(2,2)=0.0E-02_CMISSDP 804 CouplingCoeffs(2,3)=0.0E-02_CMISSDP 805 CouplingCoeffs(2,4)=0.0E-02_CMISSDP 806 CouplingCoeffs(2,5)=0.0E-02_CMISSDP 807 CouplingCoeffs(3,1)=0.0E-02_CMISSDP 808 CouplingCoeffs(3,2)=0.0E-02_CMISSDP 809 CouplingCoeffs(3,3)=0.0E-02_CMISSDP 810 CouplingCoeffs(3,4)=0.0E-02_CMISSDP 811 CouplingCoeffs(3,5)=0.0E-02_CMISSDP 812 CouplingCoeffs(4,1)=0.0E-02_CMISSDP 813 CouplingCoeffs(4,2)=0.0E-02_CMISSDP 814 CouplingCoeffs(4,3)=0.0E-02_CMISSDP 815 CouplingCoeffs(4,4)=0.0E-02_CMISSDP 816 CouplingCoeffs(4,5)=0.0E-02_CMISSDP 817 CouplingCoeffs(5,1)=0.0E-02_CMISSDP 818 CouplingCoeffs(5,2)=0.0E-02_CMISSDP 819 CouplingCoeffs(5,3)=0.0E-02_CMISSDP 820 CouplingCoeffs(5,4)=0.0E-02_CMISSDP 821 CouplingCoeffs(5,5)=0.0E-02_CMISSDP 822 ELSE 823 write(*,*) "Can't initialise coupling coefficients array." 824 ENDIF 825 !Define the material parameters for each compartments' constitutive law (for determining pressure) 826 Nparams=3 827 ALLOCATE(ConstitutiveParams(Ncompartments,Nparams)) 828 IF(Ncompartments==2)THEN 829 ConstitutiveParams(1,1)=10.0E-01_CMISSDP 830 ConstitutiveParams(1,2)=10.0E-01_CMISSDP 831 ConstitutiveParams(1,3)=10.0E-01_CMISSDP 832 ConstitutiveParams(2,1)=10.0E-01_CMISSDP 833 ConstitutiveParams(2,2)=10.0E-01_CMISSDP 834 ConstitutiveParams(2,3)=10.0E-01_CMISSDP 835 ELSE IF(Ncompartments==3)THEN 836 ConstitutiveParams(1,1)=1.0E-02_CMISSDP 837 ConstitutiveParams(1,2)=1.0E-02_CMISSDP 838 ConstitutiveParams(1,3)=0.0E-02_CMISSDP 839 ConstitutiveParams(2,1)=1.0E-02_CMISSDP 840 ConstitutiveParams(2,2)=2.0E-02_CMISSDP 841 ConstitutiveParams(2,3)=1.0E-02_CMISSDP 842 ConstitutiveParams(3,1)=0.0E-02_CMISSDP 843 ConstitutiveParams(3,2)=1.0E-02_CMISSDP 844 ConstitutiveParams(3,3)=1.0E-02_CMISSDP 845 ELSE IF(Ncompartments==4)THEN 846 ConstitutiveParams(1,1)=0.0E-02_CMISSDP 847 ConstitutiveParams(1,2)=0.0E-02_CMISSDP 848 ConstitutiveParams(1,3)=0.0E-02_CMISSDP 849 ConstitutiveParams(2,1)=0.0E-02_CMISSDP 850 ConstitutiveParams(2,2)=0.0E-02_CMISSDP 851 ConstitutiveParams(2,3)=0.0E-02_CMISSDP 852 ConstitutiveParams(3,1)=0.0E-02_CMISSDP 853 ConstitutiveParams(3,2)=0.0E-02_CMISSDP 854 ConstitutiveParams(3,3)=0.0E-02_CMISSDP 855 ConstitutiveParams(4,1)=0.0E-02_CMISSDP 856 ConstitutiveParams(4,2)=0.0E-02_CMISSDP 857 ConstitutiveParams(4,3)=0.0E-02_CMISSDP 858 ELSE IF(Ncompartments==5)THEN 859 ConstitutiveParams(1,1)=0.0E-02_CMISSDP 860 ConstitutiveParams(1,2)=0.0E-02_CMISSDP 861 ConstitutiveParams(1,3)=0.0E-02_CMISSDP 862 ConstitutiveParams(2,1)=0.0E-02_CMISSDP 863 ConstitutiveParams(2,2)=0.0E-02_CMISSDP 864 ConstitutiveParams(2,3)=0.0E-02_CMISSDP 865 ConstitutiveParams(3,1)=0.0E-02_CMISSDP 866 ConstitutiveParams(3,2)=0.0E-02_CMISSDP 867 ConstitutiveParams(3,3)=0.0E-02_CMISSDP 868 ConstitutiveParams(4,1)=0.0E-02_CMISSDP 869 ConstitutiveParams(4,2)=0.0E-02_CMISSDP 870 ConstitutiveParams(4,3)=0.0E-02_CMISSDP 871 ConstitutiveParams(5,1)=0.0E-02_CMISSDP 872 ConstitutiveParams(5,2)=0.0E-02_CMISSDP 873 ConstitutiveParams(5,3)=0.0E-02_CMISSDP 874 ELSE 875 write(*,*) "Can't initialise constitutive parameters array." 876 ENDIF 877 !MATERIALS FIELDS 878 !Auto-created field contains a U variable type to store the diffusion coefficient(s) 879 !It also contains a V variable type to store the coupling coefficients 880 DO icompartment = 1,Ncompartments 881 MaterialsFieldUserNumberDarcy = 400+icompartment 882 CALL CMISSField_Initialise(MaterialsFieldDarcy(icompartment),Err) 883 CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy(icompartment),MaterialsFieldUserNumberDarcy,& 884 & MaterialsFieldDarcy(icompartment),Err) 885 CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy(icompartment),Err) 886 END DO 887 DO icompartment = 1,Ncompartments 888 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, & 889 & CMISS_FIELD_VALUES_SET_TYPE, & 890 & MaterialsFieldUserNumberDarcyPorosity,POROSITY_PARAM_DARCY,Err) 891 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U_VARIABLE_TYPE, & 892 & CMISS_FIELD_VALUES_SET_TYPE, & 893 & MaterialsFieldUserNumberDarcyPermOverVis,PERM_OVER_VIS_PARAM_DARCY,Err) 894 END DO 895 DO icompartment = 1, Ncompartments 896 DO COMPONENT_NUMBER=1, Ncompartments 897 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, & 898 & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err) 899! CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, & 900! & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err) 901 END DO 902 END DO 903 DO icompartment = 1, Ncompartments 904 DO COMPONENT_NUMBER=1,Nparams 905 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy(icompartment),CMISS_FIELD_U1_VARIABLE_TYPE, & 906 & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,ConstitutiveParams(icompartment,COMPONENT_NUMBER),Err) 907! CALL CMISSField_ParameterSetUpdateConstant(MaterialsFieldDarcy(icompartment),CMISS_FIELD_V_VARIABLE_TYPE, & 908! & CMISS_FIELD_VALUES_SET_TYPE,COMPONENT_NUMBER,CouplingCoeffs(icompartment,COMPONENT_NUMBER),Err) 909 END DO 910 END DO 911 912 ! 913 !================================================================================================================================ 914 ! 915 916 !EQUATIONS SET EQUATIONS 917 918 !Darcy 919 DO icompartment=1,Ncompartments 920 !Create the equations set equations 921 CALL CMISSEquations_Initialise(EquationsDarcy(icompartment),Err) 922 CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy(icompartment),EquationsDarcy(icompartment),Err) 923 !Set the equations matrices sparsity type 924 CALL CMISSEquations_SparsityTypeSet(EquationsDarcy(icompartment…
Large files files are truncated, but you can click here to view the full file