/MultiPhysics/Poroelasticity/FiniteElasticityDarcy/QuadraticEllipsoidDrivenDarcy/src/QuadraticEllipsoidDrivenDarcyExample.f90
FORTRAN Modern | 1124 lines | 649 code | 181 blank | 294 comment | 0 complexity | 065c8a8617ed85fdc58477d6f0c89bff 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 QUADRATICELLIPSOIDDRIVENDARCYEXAMPLE 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=2 ! 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 :: FieldDependentUserNumber=4 134 INTEGER(CMISSIntg), PARAMETER :: FieldDependentNumberOfVariables=4 !2 135 INTEGER(CMISSIntg), PARAMETER :: FieldDependentSolidNumberOfComponents=4 136 INTEGER(CMISSIntg), PARAMETER :: FieldDependentFluidNumberOfComponents=4 137 138 INTEGER(CMISSIntg), PARAMETER :: IndependentFieldDarcyUserNumber=15 139 140 INTEGER(CMISSIntg), PARAMETER :: EquationSetUserNumberSolid=1 141 INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumberSolid=13 142 INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=1 143 144 145 INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberDarcy=8 146 INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberDarcy=12 147 INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumberDarcy=22 148 INTEGER(CMISSIntg), PARAMETER :: SourceFieldDarcyUserNumber=42 149 150 INTEGER(CMISSIntg), PARAMETER :: ControlLoopSolidNumber=1 151 INTEGER(CMISSIntg), PARAMETER :: ControlLoopFluidNumber=2 152 INTEGER(CMISSIntg), PARAMETER :: ControlLoopSubiterationNumber=1 153 INTEGER(CMISSIntg), PARAMETER :: SolverSolidIndex=1 154 INTEGER(CMISSIntg), PARAMETER :: SolverDarcyIndex=1 155 156 !Program types 157 158 159 !Program variables 160 161 INTEGER(CMISSIntg) :: MPI_IERROR 162 INTEGER(CMISSIntg) :: EquationsSetIndex 163 INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber 164 REAL(CMISSDP) :: FibreFieldAngle(3) 165 REAL(CMISSDP) :: nu,theta,omega,XI3,XI3delta,XI2delta, zero 166 INTEGER(CMISSIntg) ::i,j,k,component_idx,node_idx,TOTAL_NUMBER_NODES_XI(3) 167 !For grabbing surfaces 168 INTEGER(CMISSIntg) :: InnerNormalXi,OuterNormalXi,TopNormalXi 169 INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodes(:) 170 INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodes(:) 171 INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodes(:) 172 INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodesDarcyVel(:) 173 INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodesDarcyVel(:) 174 INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodesDarcyVel(:) 175 INTEGER(CMISSIntg) :: NN,NODE,NodeDomain 176 REAL(CMISSDP) :: XCoord,YCoord,ZCoord 177 LOGICAL :: X_FIXED,Y_FIXED,X_OKAY,Y_OKAY 178 179 INTEGER(CMISSIntg) :: GeometricFieldDarcyMeshComponentNumber, DarcyVelMeshComponentNumber, DarcyMassIncreaseMeshComponentNumber 180 INTEGER(CMISSIntg) :: MeshComponentNumber_dummy 181 182 INTEGER(CMISSIntg) :: NUMBER_OF_DOMAINS 183 184 INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS 185 186 INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS 187 INTEGER(CMISSIntg) :: RESTART_VALUE 188 189 INTEGER(CMISSIntg) :: EQUATIONS_DARCY_OUTPUT 190 INTEGER(CMISSIntg) :: COMPONENT_NUMBER, NODE_NUMBER 191 INTEGER(CMISSIntg) :: CONDITION 192 193 INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY 194 INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE 195 INTEGER(CMISSIntg) :: LINEAR_SOLVER_DARCY_OUTPUT_TYPE 196 197 REAL(CMISSDP) :: COORD_X, COORD_Y, COORD_Z 198 REAL(CMISSDP) :: DOMAIN_X1, DOMAIN_X2, DOMAIN_Y1, DOMAIN_Y2, DOMAIN_Z1, DOMAIN_Z2 199 REAL(CMISSDP) :: GEOMETRY_TOLERANCE 200 INTEGER(CMISSIntg) :: EDGE_COUNT 201 INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SOLID 202 REAL(CMISSDP) :: INITIAL_FIELD_DARCY(4) 203 REAL(CMISSDP) :: DIVERGENCE_TOLERANCE 204 REAL(CMISSDP) :: RELATIVE_TOLERANCE 205 REAL(CMISSDP) :: ABSOLUTE_TOLERANCE 206 REAL(CMISSDP) :: LINESEARCH_ALPHA 207 REAL(CMISSDP) :: VALUE 208 REAL(CMISSDP) :: POROSITY_PARAM_DARCY, PERM_OVER_VIS_PARAM_DARCY 209 210 LOGICAL :: EXPORT_FIELD_IO 211 LOGICAL :: LINEAR_SOLVER_DARCY_DIRECT_FLAG 212 213 !CMISS variables 214 215 TYPE(CMISSBasisType) :: QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis 216 TYPE(CMISSBoundaryConditionsType) :: BoundaryConditions 217 TYPE(CMISSCoordinateSystemType) :: CoordinateSystem, WorldCoordinateSystem 218 TYPE(CMISSMeshType) :: Mesh 219 TYPE(CMISSGeneratedMeshType) :: GeneratedMesh 220 TYPE(CMISSDecompositionType) :: Decomposition 221 TYPE(CMISSEquationsType) :: Equations 222 TYPE(CMISSEquationsSetType) :: EquationsSetSolid 223 TYPE(CMISSFieldType) :: GeometricFieldSolid,FibreFieldSolid,MaterialFieldSolid 224 TYPE(CMISSFieldType) :: DependentField,EquationsSetFieldSolid 225 226 TYPE(CMISSFieldType) :: IndependentFieldDarcy 227 228 TYPE(CMISSFieldType) :: GeometricFieldDarcy 229 230 TYPE(CMISSFieldType) :: MaterialsFieldDarcy 231 TYPE(CMISSFieldType) :: EquationsSetFieldDarcy 232 TYPE(CMISSFieldType) :: SourceFieldDarcy 233 !Boundary conditions 234 TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsDarcy 235 !Equations sets 236 TYPE(CMISSEquationsSetType) :: EquationsSetDarcy 237 !Equations 238 TYPE(CMISSEquationsType) :: EquationsDarcy 239 240 TYPE(CMISSFieldsType) :: Fields 241 TYPE(CMISSProblemType) :: Problem 242 TYPE(CMISSRegionType) :: Region,WorldRegion 243 TYPE(CMISSSolverType) :: SolverSolid,LinearSolverSolid 244 TYPE(CMISSSolverEquationsType) :: SolverEquationsSolid 245 TYPE(CMISSControlLoopType) :: ControlLoop 246 247 !Solvers 248 TYPE(CMISSSolverType) :: DynamicSolverDarcy 249 TYPE(CMISSSolverType) :: LinearSolverDarcy 250 !Solver equations 251 TYPE(CMISSSolverEquationsType) :: SolverEquationsDarcy 252 253 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_START_TIME 254 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_STOP_TIME 255 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_THETA 256 REAL(CMISSDP) :: DYNAMIC_SOLVER_DARCY_TIME_INCREMENT 257 258 INTEGER(CMISSIntg),allocatable :: ElementUserNodes(:) 259 INTEGER(CMISSIntg) :: NUMBER_USER_ELEMENT_NODES, ELEMENT_NUMBER 260 261#ifdef WIN32 262 !Quickwin type 263 LOGICAL :: QUICKWIN_STATUS=.FALSE. 264 TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG 265#endif 266 267 !Generic CMISS variables 268 INTEGER(CMISSIntg) :: Err 269 270#ifdef WIN32 271 !Initialise QuickWin 272 QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title 273 QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows 274 QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN 275 !Set the window parameters 276 QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG) 277 !If attempt fails set with system estimated values 278 IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG) 279#endif 280 281 !Set initial values 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=1.0E-3_CMISSDP 298 DYNAMIC_SOLVER_DARCY_TIME_INCREMENT=1.0e-3_CMISSDP 299 DYNAMIC_SOLVER_DARCY_STOP_TIME=2_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 312 313 !LinearMeshComponentNumber/QuadraticMeshComponentNumber 314 DarcyVelMeshComponentNumber = LinearMeshComponentNumber 315 DarcyMassIncreaseMeshComponentNumber = LinearMeshComponentNumber 316! GeometricFieldDarcyMeshComponentNumber = DarcyVelMeshComponentNumber 317 GeometricFieldDarcyMeshComponentNumber = QuadraticMeshComponentNumber 318 319 320 ! 321 !================================================================================================================================ 322 ! 323 324 !Intialise cmiss 325 CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err) 326 327 CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err) 328 329 WRITE(*,'(A)') "Program starting." 330 331 !Set all diganostic levels on for testing 332 CALL CMISSDiagnosticsSetOn(CMISS_FROM_DIAG_TYPE,(/1,2,3,4,5/),"Diagnostics",(/"PROBLEM_FINITE_ELEMENT_CALCULATE"/),Err) 333 334 !Get the number of computational nodes and this computational node number 335 CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err) 336 CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err) 337 338 !Broadcast the number of elements in the X,Y and Z directions and the number of partitions to the other computational nodes 339! CALL MPI_BCAST(NumberGlobalXElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 340! CALL MPI_BCAST(NumberGlobalYElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 341! CALL MPI_BCAST(NumberGlobalZElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 342! CALL MPI_BCAST(NumberOfDomains,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR) 343 NumberOfDomains=NumberOfComputationalNodes 344 345 !Create a CS - default is 3D rectangular cartesian CS with 0,0,0 as origin 346 CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err) 347 CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err) 348 CALL CMISSCoordinateSystem_TypeSet(CoordinateSystem,CMISS_COORDINATE_RECTANGULAR_CARTESIAN_TYPE,Err) 349 CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NumberOfSpatialCoordinates,Err) 350 CALL CMISSCoordinateSystem_OriginSet(CoordinateSystem,(/0.0_CMISSDP,0.0_CMISSDP,0.0_CMISSDP/),Err) 351 CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err) 352 353 ! 354 !================================================================================================================================ 355 ! 356 357 !Create a region and assign the CS to the region 358 CALL CMISSRegion_Initialise(Region,Err) 359 CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err) 360 CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err) 361 CALL CMISSRegion_CreateFinish(Region,Err) 362 363 ! 364 !================================================================================================================================ 365 ! 366 367 !Define basis functions - tri-linear Lagrange and tri-Quadratic Lagrange, each with collapsed variant 368 !Quadratic Basis 369 CALL CMISSBasis_Initialise(QuadraticBasis,Err) 370 CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err) 371 CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, & 372 & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err) 373 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, & 374 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 375! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 376 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err) !Have to do this 377 CALL CMISSBasis_CreateFinish(QuadraticBasis,Err) 378 379 !Collapsed Quadratic Basis 380 CALL CMISSBasis_Initialise(QuadraticCollapsedBasis,Err) 381 CALL CMISSBasis_CreateStart(QuadraticCollapsedBasisUserNumber,QuadraticCollapsedBasis,Err) 382 CALL CMISSBasis_TypeSet(QuadraticCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err) 383 CALL CMISSBasis_NumberOfXiSet(QuadraticCollapsedBasis,NumberOfXiCoordinates,Err) 384 CALL CMISSBasis_InterpolationXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, & 385 & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err) 386 CALL CMISSBasis_CollapsedXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED, & 387 & CMISS_BASIS_COLLAPSED_AT_XI0,CMISS_BASIS_NOT_COLLAPSED/),Err) 388 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticCollapsedBasis, & 389 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 390! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 391 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticCollapsedBasis,.true.,Err) !Have to do this 392 CALL CMISSBasis_CreateFinish(QuadraticCollapsedBasis,Err) 393 394 !Linear Basis 395 CALL CMISSBasis_Initialise(LinearBasis,Err) 396 CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err) 397 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, & 398 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 399! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 400 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err) !Have to do this (unused) due to field_interp setup 401 CALL CMISSBasis_CreateFinish(LinearBasis,Err) 402 403 !Collapsed Linear Basis 404 CALL CMISSBasis_Initialise(LinearCollapsedBasis,Err) 405 CALL CMISSBasis_CreateStart(LinearCollapsedBasisUserNumber,LinearCollapsedBasis,Err) 406 CALL CMISSBasis_TypeSet(LinearCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err) 407 CALL CMISSBasis_NumberOfXiSet(LinearCollapsedBasis,NumberOfXiCoordinates,Err) 408 CALL CMISSBasis_InterpolationXiSet(LinearCollapsedBasis,(/CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION, & 409 & CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION,CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION/),Err) 410 CALL CMISSBasis_CollapsedXiSet(LinearCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED,CMISS_BASIS_COLLAPSED_AT_XI0, & 411 & CMISS_BASIS_NOT_COLLAPSED/),Err) 412 CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearCollapsedBasis, & 413 & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err) 414! & (/CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME,CMISS_BASIS_HIGH_QUADRATURE_SCHEME/),Err) 415 CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearCollapsedBasis,.true.,Err) !Have to do this (unused) due to field_interp setup 416 CALL CMISSBasis_CreateFinish(LinearCollapsedBasis,Err) 417 418 ! 419 !================================================================================================================================ 420 ! 421 422 !Start the creation of a generated ellipsoid mesh 423 CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err) 424 CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err) 425 !Set up an ellipsoid mesh 426 CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_MESH_TYPE,Err) 427 !Set the quadratic and linear bases 428 CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,[QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis],Err) 429 !Define the mesh on the region 430 CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/LONG_AXIS,SHORT_AXIS,WALL_THICKNESS,CUTOFF_ANGLE/),Err) 431 CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NumberGlobalXElements,NumberGlobalYElements, & 432 & NumberGlobalZElements/),Err) 433 434 !Finish the creation of a generated mesh in the region 435 CALL CMISSMesh_Initialise(Mesh,Err) 436 CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err) 437 438 ! 439 !================================================================================================================================ 440 ! 441 442 !Create a decomposition 443 CALL CMISSDecomposition_Initialise(Decomposition,Err) 444 CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err) 445 CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err) 446 CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err) 447 CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.TRUE.,Err) 448 CALL CMISSDecomposition_CreateFinish(Decomposition,Err) 449 450 ! 451 !================================================================================================================================ 452 ! 453 454 ! --- GeometricFieldSolid --- 455 !Create a field to put the geometry (default is geometry) 456 CALL CMISSField_Initialise(GeometricFieldSolid,Err) 457 CALL CMISSField_CreateStart(FieldGeometryUserNumberSolid,Region,GeometricFieldSolid,Err) 458 CALL CMISSField_MeshDecompositionSet(GeometricFieldSolid,Decomposition,Err) 459 CALL CMISSField_TypeSet(GeometricFieldSolid,CMISS_FIELD_GEOMETRIC_TYPE,Err) 460 CALL CMISSField_NumberOfVariablesSet(GeometricFieldSolid,FieldGeometryNumberOfVariables,Err) 461 CALL CMISSField_NumberOfComponentsSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err) 462 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err) 463 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err) 464 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err) 465 CALL CMISSField_CreateFinish(GeometricFieldSolid,Err) 466 467 !Update the geometric field parameters 468 CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldSolid,Err) 469 470 ! 471 !================================================================================================================================ 472 ! 473 474 ! --- GeometricFieldDarcy --- 475 !Create a field to put the geometry (default is geometry) 476 CALL CMISSField_Initialise(GeometricFieldDarcy,Err) 477 CALL CMISSField_CreateStart(FieldGeometryUserNumberDarcy,Region,GeometricFieldDarcy,Err) 478 CALL CMISSField_MeshDecompositionSet(GeometricFieldDarcy,Decomposition,Err) 479 CALL CMISSField_TypeSet(GeometricFieldDarcy,CMISS_FIELD_GEOMETRIC_TYPE,Err) 480 CALL CMISSField_NumberOfVariablesSet(GeometricFieldDarcy,FieldGeometryNumberOfVariables,Err) 481 CALL CMISSField_NumberOfComponentsSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err) 482 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,1, & 483 & GeometricFieldDarcyMeshComponentNumber,Err) 484 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,2, & 485 & GeometricFieldDarcyMeshComponentNumber,Err) 486 CALL CMISSField_ComponentMeshComponentSet(GeometricFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,3, & 487 & GeometricFieldDarcyMeshComponentNumber,Err) 488 CALL CMISSField_CreateFinish(GeometricFieldDarcy,Err) 489 490 !Update the geometric field parameters 491 CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricFieldDarcy,Err) 492 493 ! 494 !================================================================================================================================ 495 ! 496 497 !Create a fibre field and attach it to the geometric field 498 CALL CMISSField_Initialise(FibreFieldSolid,Err) 499 CALL CMISSField_CreateStart(FieldFibreUserNumber,Region,FibreFieldSolid,Err) 500 CALL CMISSField_TypeSet(FibreFieldSolid,CMISS_FIELD_FIBRE_TYPE,Err) 501 CALL CMISSField_MeshDecompositionSet(FibreFieldSolid,Decomposition,Err) 502 CALL CMISSField_GeometricFieldSet(FibreFieldSolid,GeometricFieldSolid,Err) 503 CALL CMISSField_NumberOfVariablesSet(FibreFieldSolid,FieldFibreNumberOfVariables,Err) 504 CALL CMISSField_NumberOfComponentsSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreNumberOfComponents,Err) 505 CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err) 506 CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err) 507 CALL CMISSField_ComponentMeshComponentSet(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err) 508 CALL CMISSField_CreateFinish(FibreFieldSolid,Err) 509 510 !Set Fibre directions (this block is parallel-untested) 511 node_idx=0 512 !This is valid only for quadratic basis functions 513 TOTAL_NUMBER_NODES_XI(1)=NumberGlobalXElements*2 514 TOTAL_NUMBER_NODES_XI(2)=NumberGlobalYElements*2+1 515 TOTAL_NUMBER_NODES_XI(3)=NumberGlobalZElements*2+1 516 517 XI2delta=(PI-CUTOFF_ANGLE)/(TOTAL_NUMBER_NODES_XI(2)-1) 518 XI3=0 519 XI3delta=(1.0)/(TOTAL_NUMBER_NODES_XI(3)-1) 520 zero=0 521 DO k=1, TOTAL_NUMBER_NODES_XI(3) 522 !Apex nodes 523 j=1 524 i=1 525 node_idx=node_idx+1 526 CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err) 527 IF(NodeDomain==ComputationalNodeNumber) THEN 528 FibreFieldAngle=(/zero,zero,zero/) 529 DO component_idx=1,FieldFibreNumberOfComponents 530 CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, & 531 & DerivativeUserNumber,node_idx,component_idx,FibreFieldAngle(component_idx),Err) 532 ENDDO 533 ENDIF 534 theta=atan(FIBRE_SLOPE_CHANGE*XI3+FIBRE_SLOPE_INTERSECTION) 535 DO j=2, TOTAL_NUMBER_NODES_XI(2) 536 nu=PI-XI2delta*(j-1) 537 omega=PI/2+cos(2*nu)*atan(SHEET_SLOPE_BASE_ENDO)*(-2*XI3+1) 538 DO i=1, TOTAL_NUMBER_NODES_XI(1) 539 node_idx=node_idx+1 540 CALL CMISSDecomposition_NodeDomainGet(Decomposition,node_idx,1,NodeDomain,Err) 541 IF(NodeDomain==ComputationalNodeNumber) THEN 542 FibreFieldAngle=(/theta,zero,omega/) 543 DO component_idx=1,FieldFibreNumberOfComponents 544 CALL CMISSField_ParameterSetUpdateNode(FibreFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, & 545 & DerivativeUserNumber, node_idx,component_idx,FibreFieldAngle(component_idx),Err) 546 ENDDO 547 ENDIF 548 ENDDO 549 ENDDO 550 XI3=XI3+XI3delta 551 ENDDO 552 553 !Create a material field and attach it to the geometric field 554 CALL CMISSField_Initialise(MaterialFieldSolid,Err) 555 CALL CMISSField_CreateStart(FieldMaterialUserNumber,Region,MaterialFieldSolid,Err) 556 CALL CMISSField_TypeSet(MaterialFieldSolid,CMISS_FIELD_MATERIAL_TYPE,Err) 557 CALL CMISSField_MeshDecompositionSet(MaterialFieldSolid,Decomposition,Err) 558 CALL CMISSField_GeometricFieldSet(MaterialFieldSolid,GeometricFieldSolid,Err) 559 CALL CMISSField_NumberOfVariablesSet(MaterialFieldSolid,FieldMaterialNumberOfVariables,Err) 560 CALL CMISSField_NumberOfComponentsSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialNumberOfComponents,Err) 561 CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_CONSTANT_INTERPOLATION,Err) 562 CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_CONSTANT_INTERPOLATION,Err) 563 CALL CMISSField_ComponentInterpolationSet(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,3,CMISS_FIELD_CONSTANT_INTERPOLATION,Err) 564 CALL CMISSField_CreateFinish(MaterialFieldSolid,Err) 565 566 !Set Mooney-Rivlin constants 567 CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,C1,Err) 568 CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,C2,Err) 569 CALL CMISSField_ComponentValuesInitialise(MaterialFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,C3,Err) 570 571 ! 572 !================================================================================================================================ 573 ! 574 575 !EQUATIONS SETS 576 577 !Create the equations set for ALE Darcy 578 CALL CMISSField_Initialise(EquationsSetFieldDarcy,Err) 579 CALL CMISSEquationsSet_Initialise(EquationsSetDarcy,Err) 580 CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberDarcy,Region,GeometricFieldDarcy, & 581 & CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, & 582 & CMISS_EQUATIONS_SET_DARCY_EQUATION_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE,& 583 & EquationsSetFieldUserNumberDarcy,EquationsSetFieldDarcy,EquationsSetDarcy,Err) 584 CALL CMISSEquationsSet_CreateFinish(EquationsSetDarcy,Err) 585 586 !Create the equations set for the solid 587 CALL CMISSField_Initialise(EquationsSetFieldSolid,Err) 588 CALL CMISSEquationsSet_Initialise(EquationsSetSolid,Err) 589 CALL CMISSEquationsSet_CreateStart(EquationSetUserNumberSolid,Region,FibreFieldSolid,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, & 590 & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_INCOMPRESS_ELASTICITY_DRIVEN_DARCY_SUBTYPE, & 591 & EquationsSetFieldUserNumberSolid,EquationsSetFieldSolid,EquationsSetSolid,Err) 592 CALL CMISSEquationsSet_CreateFinish(EquationsSetSolid,Err) 593 594 ! 595 !================================================================================================================================ 596 ! 597 598 !DEPENDENT FIELDS 599 600 !Create a dependent field with four variables (U, DelUDelN = solid, V, DelVDelN = Darcy) and four components 601 !Solid: The U, DelUDelN variables have 4 components (3 displacement, 1 pressure) 602 CALL CMISSField_Initialise(DependentField,Err) 603 CALL CMISSField_CreateStart(FieldDependentUserNumber,Region,DependentField,Err) 604 CALL CMISSField_TypeSet(DependentField,CMISS_FIELD_GENERAL_TYPE,Err) 605 CALL CMISSField_MeshDecompositionSet(DependentField,Decomposition,Err) 606 CALL CMISSField_GeometricFieldSet(DependentField,GeometricFieldSolid,Err) 607 CALL CMISSField_DependentTypeSet(DependentField,CMISS_FIELD_DEPENDENT_TYPE,Err) 608 CALL CMISSField_NumberOfVariablesSet(DependentField,FieldDependentNumberOfVariables,Err) 609 610 CALL CMISSField_VariableTypesSet(DependentField,(/CMISS_FIELD_U_VARIABLE_TYPE, & 611 & CMISS_FIELD_DELUDELN_VARIABLE_TYPE,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_DELVDELN_VARIABLE_TYPE/),Err) 612 CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err) 613 CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,FieldDependentSolidNumberOfComponents,Err) 614 CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err) 615 CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,FieldDependentFluidNumberOfComponents,Err) 616 617 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err) 618 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err) 619 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err) 620 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err) 621 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err) 622 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err) 623 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err) 624 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err) 625 626 !Darcy: The V, DelVDelN variables have 4 components (3 velocities, 1 mass increase) 627 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber,Err) 628 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber,Err) 629 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber,Err) 630 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,4,DarcyMassIncreaseMeshComponentNumber,Err) 631 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber,Err) 632 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber,Err) 633 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber,Err) 634 CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELVDELN_VARIABLE_TYPE,4, & 635 & DarcyMassIncreaseMeshComponentNumber,Err) 636 637 CALL CMISSField_ScalingTypeSet(DependentField,CMISS_FIELD_UNIT_SCALING,Err) 638 639 CALL CMISSField_CreateFinish(DependentField,Err) 640 641 ! 642 !================================================================================================================================ 643 ! 644 645 CALL CMISSEquationsSet_DependentCreateStart(EquationsSetSolid,FieldDependentUserNumber,DependentField,Err) 646 CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetSolid,Err) 647 648 CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetSolid,FieldMaterialUserNumber,MaterialFieldSolid,Err) 649 CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetSolid,Err) 650 651 ! 652 !================================================================================================================================ 653 ! 654 655 CALL CMISSEquationsSet_DependentCreateStart(EquationsSetDarcy,FieldDependentUserNumber,DependentField,Err) 656 CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetDarcy,Err) 657 658 DO COMPONENT_NUMBER=1,FieldDependentFluidNumberOfComponents 659 CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_V_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 660 & COMPONENT_NUMBER,INITIAL_FIELD_DARCY(COMPONENT_NUMBER),Err) 661 ENDDO 662 663 ! 664 !================================================================================================================================ 665 ! 666 667 !INDEPENDENT FIELD Darcy for storing BC flags 668 669 CALL CMISSField_Initialise(IndependentFieldDarcy,Err) 670 CALL CMISSEquationsSet_IndependentCreateStart(EquationsSetDarcy,IndependentFieldDarcyUserNumber, & 671 & IndependentFieldDarcy,Err) 672 673 CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,1,DarcyVelMeshComponentNumber,Err) 674 CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,2,DarcyVelMeshComponentNumber,Err) 675 CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,3,DarcyVelMeshComponentNumber,Err) 676! CALL CMISSField_ComponentMeshComponentSet(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,4,DarcyMassIncreaseMeshComponentNumber,Err) 677 678 CALL CMISSEquationsSet_IndependentCreateFinish(EquationsSetDarcy,Err) 679 680 CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, & 681 & 0.0_CMISSDP,Err) 682 CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2, & 683 & 0.0_CMISSDP,Err) 684 CALL CMISSField_ComponentValuesInitialise(IndependentFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3, & 685 & 0.0_CMISSDP,Err) 686 687 ! 688 !================================================================================================================================ 689 ! 690 691 !Create the equations set materials field variables for ALE Darcy 692 CALL CMISSField_Initialise(MaterialsFieldDarcy,Err) 693 CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetDarcy,MaterialsFieldUserNumberDarcy, & 694 & MaterialsFieldDarcy,Err) 695 CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetDarcy,Err) 696 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 697 & 1,POROSITY_PARAM_DARCY,Err) 698 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 699 & 2,PERM_OVER_VIS_PARAM_DARCY,Err) 700 701 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 702 & 3,0.0_CMISSDP,Err) 703 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 704 & 4,0.0_CMISSDP,Err) 705 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 706 & 5,PERM_OVER_VIS_PARAM_DARCY,Err) 707 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 708 & 6,0.0_CMISSDP,Err) 709 CALL CMISSField_ComponentValuesInitialise(MaterialsFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 710 & 7,PERM_OVER_VIS_PARAM_DARCY,Err) 711 712 ! 713 !================================================================================================================================ 714 ! 715 716 !Source field 717 CALL CMISSField_Initialise(SourceFieldDarcy,Err) 718 CALL CMISSEquationsSet_SourceCreateStart(EquationsSetDarcy,SourceFieldDarcyUserNumber,SourceFieldDarcy,Err) 719 CALL CMISSEquationsSet_SourceCreateFinish(EquationsSetDarcy,Err) 720 721! ELEMENT_NUMBER = 5 722! COMPONENT_NUMBER = 4 723! VALUE = 4.2_CMISSDP 724! ! CALL CMISSField_ParameterSetUpdateElement(RegionUserNumber,SourceFieldDarcyUserNumber,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 725! ! & ELEMENT_NUMBER,COMPONENT_NUMBER,VALUE,Err) 726! 727! NUMBER_USER_ELEMENT_NODES = 27 !hardcoding is bad - but how to access the number of element nodes ? 728! !- there is no CMISS library call available for this 729! !- traversing the structure of 'CMISSMeshElementsType' does not work either, 730! ! since certain members are private 731! 732! allocate( ElementUserNodes(NUMBER_USER_ELEMENT_NODES) ) 733! 734! CALL CMISSMeshElements_NodesGet(RegionUserNumber,MeshUserNumber,QuadraticMeshComponentNumber,ELEMENT_NUMBER, & 735! & ElementUserNodes,Err) 736! 737! DO NN=1,NUMBER_USER_ELEMENT_NODES 738! NODE_NUMBER = ElementUserNodes(NN) 739! CALL CMISSField_ParameterSetUpdateNode(SourceFieldDarcy,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 740! & 1,NODE_NUMBER,COMPONENT_NUMBER,VALUE,Err) 741! ENDDO 742 743 ! 744 !================================================================================================================================ 745 ! 746 747 !EQUATIONS SET EQUATIONS 748 749 !Darcy 750 CALL CMISSEquations_Initialise(EquationsDarcy,Err) 751 CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetDarcy,EquationsDarcy,Err) 752 CALL CMISSEquations_SparsityTypeSet(EquationsDarcy,CMISS_EQUATIONS_SPARSE_MATRICES,Err) 753 CALL CMISSEquations_OutputTypeSet(EquationsDarcy,EQUATIONS_DARCY_OUTPUT,Err) 754 CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetDarcy,Err) 755 756 !Solid 757 CALL CMISSEquations_Initialise(Equations,Err) 758 CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetSolid,Equations,Err) 759 CALL CMISSEquations_SparsityTypeSet(Equations,CMISS_EQUATIONS_SPARSE_MATRICES,Err) 760 CALL CMISSEquations_OutputTypeSet(Equations,CMISS_EQUATIONS_NO_OUTPUT,Err) 761 CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetSolid,Err) 762 763 ! 764 !================================================================================================================================ 765 ! 766 767 !Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure 768 CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, & 769 & CMISS_FIELD_VALUES_SET_TYPE, & 770 & 1,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,Err) 771 CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, & 772 & CMISS_FIELD_VALUES_SET_TYPE, & 773 & 2,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,Err) 774 CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE, & 775 & CMISS_FIELD_VALUES_SET_TYPE, & 776 & 3,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,Err) 777! CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4,-14.0_CMISSDP,Err) 778 CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4,0.0_CMISSDP, & 779 & Err) 780 781 ! 782 !================================================================================================================================ 783 ! 784 785 !Define the problem 786 CALL CMISSProblem_Initialise(Problem,Err) 787 CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err) 788 CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_MULTI_PHYSICS_CLASS,CMISS_PROBLEM_FINITE_ELASTICITY_DARCY_TYPE, & 789 & CMISS_PROBLEM_QUASISTATIC_ELASTICITY_TRANSIENT_DARCY_SUBTYPE,Err) 790 CALL CMISSProblem_CreateFinish(Problem,Err) 791 792 ! 793 !================================================================================================================================ 794 ! 795 796 !Create the problem control loop 797 CALL CMISSProblem_ControlLoopCreateStart(Problem,Err) 798 CALL CMISSControlLoop_Initialise(ControlLoop,Err) 799 CALL CMISSProblem_ControlLoopGet(Problem,CMISS_CONTROL_LOOP_NODE,ControlLoop,Err) 800! CALL CMISSControlLoop_MaximumIterationsSet(ControlLoop,1,Err) ! this one sets the increment loop counter 801 CALL CMISSControlLoop_TimesSet(ControlLoop,DYNAMIC_SOLVER_DARCY_START_TIME,DYNAMIC_SOLVER_DARCY_STOP_TIME, & 802 & DYNAMIC_SOLVER_DARCY_TIME_INCREMENT,Err) 803 CALL CMISSControlLoop_TimeOutputSet(ControlLoop,DYNAMIC_SOLVER_DARCY_OUTPUT_FREQUENCY,Err) 804! CALL CMISSControlLoop_OutputTypeSet(ControlLoop,CMISS_CONTROL_LOOP_PROGRESS_OUTPUT,Err) 805 CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err) 806 807 ! 808 !================================================================================================================================ 809 ! 810 811 !Create the problem solvers 812 CALL CMISSSolver_Initialise(SolverSolid,Err) 813 CALL CMISSSolver_Initialise(LinearSolverSolid,Err) 814 CALL CMISSSolver_Initialise(DynamicSolverDarcy,Err) 815 CALL CMISSSolver_Initialise(LinearSolverDarcy,Err) 816 817 CALL CMISSProblem_SolversCreateStart(Problem,Err) 818 819 ! Solid 820 CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), & 821 & SolverSolidIndex,SolverSolid,Err) 822 CALL CMISSSolver_OutputTypeSet(SolverSolid,CMISS_SOLVER_PROGRESS_OUTPUT,Err) 823! CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_FD_CALCULATED,Err) 824 CALL CMISSSolver_NewtonJacobianCalculationTypeSet(SolverSolid,CMISS_SOLVER_NEWTON_JACOBIAN_EQUATIONS_CALCULATED,Err) 825 826! CALL CMISSSolverNonLinearTypeSet(SolverSolid,CMISS_SOLVER_NONLINEAR_NEWTON,Err) 827! CALL CMISSSolver_LibraryTypeSet(SolverSolid,CMISS_SOLVER_PETSC_LIBRARY,Err) 828 829 CALL CMISSSolver_NewtonAbsoluteToleranceSet(SolverSolid,ABSOLUTE_TOLERANCE,Err) 830 CALL CMISSSolver_NewtonRelativeToleranceSet(SolverSolid,RELATIVE_TOLERANCE,Err) 831 CALL CMISSSolver_NewtonMaximumIterationsSet(SolverSolid,MAXIMUM_ITERATIONS,Err) 832 833 CALL CMISSSolver_NewtonLinearSolverGet(SolverSolid,LinearSolverSolid,Err) 834 CALL CMISSSolver_LinearTypeSet(LinearSolverSolid,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err) 835 836 !Darcy 837 CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), & 838 & SolverDarcyIndex,DynamicSolverDarcy,Err) 839 CALL CMISSSolver_OutputTypeSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_OUTPUT_TYPE,Err) 840 CALL CMISSSolver_DynamicThetaSet(DynamicSolverDarcy,DYNAMIC_SOLVER_DARCY_THETA,Err) 841! CALL CMISSSolverDynamicDynamicSet(DynamicSolverDarcy,.TRUE.,Err) 842 CALL CMISSSolver_DynamicLinearSolverGet(DynamicSolverDarcy,LinearSolverDarcy,Err) 843 IF(LINEAR_SOLVER_DARCY_DIRECT_FLAG) THEN 844 CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err) 845 CALL CMISSSolver_LibraryTypeSet(LinearSolverDarcy,CMISS_SOLVER_MUMPS_LIBRARY,Err) 846 ELSE 847 CALL CMISSSolver_LinearTypeSet(LinearSolverDarcy,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err) 848 CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverDarcy,MAXIMUM_ITERATIONS,Err) 849 CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverDarcy,DIVERGENCE_TOLERANCE,Err) 850 CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverDarcy,RELATIVE_TOLERANCE,Err) 851 CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverDarcy,ABSOLUTE_TOLERANCE,Err) 852 CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverDarcy,RESTART_VALUE,Err) 853 ENDIF 854 855CALL CMISSProblem_SolversCreateFinish(Problem,Err) 856 857 858 ! 859 !================================================================================================================================ 860 ! 861 862 !Create the problem solver equations 863 CALL CMISSSolver_Initialise(SolverSolid,Err) 864 CALL CMISSSolver_Initialise(LinearSolverDarcy,Err) 865 866 CALL CMISSSolverEquations_Initialise(SolverEquationsSolid,Err) 867 CALL CMISSSolverEquations_Initialise(SolverEquationsDarcy,Err) 868 869 CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err) 870 ! 871 !Get the finite elasticity solver equations 872 CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopSolidNumber,CMISS_CONTROL_LOOP_NODE/), & 873 & SolverSolidIndex,SolverSolid,Err) 874 CALL CMISSSolver_SolverEquationsGet(SolverSolid,SolverEquationsSolid,Err) 875 CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsSolid,CMISS_SOLVER_SPARSE_MATRICES,Err) 876 CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsSolid,EquationsSetSolid,EquationsSetIndex,Err) 877 ! 878 !Get the Darcy solver equations 879 CALL CMISSProblem_SolverGet(Problem,(/ControlLoopSubiterationNumber,ControlLoopFluidNumber,CMISS_CONTROL_LOOP_NODE/), & 880 & SolverDarcyIndex,LinearSolverDarcy,Err) 881 CALL CMISSSolver_SolverEquationsGet(LinearSolverDarcy,SolverEquationsDarcy,Err) 882 CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsDarcy,CMISS_SOLVER_SPARSE_MATRICES,Err) 883 CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsDarcy,EquationsSetDarcy,EquationsSetIndex,Err) 884 ! 885 CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err) 886 887 ! 888 !================================================================================================================================ 889 ! 890 891 !Prescribe boundary conditions (absolute nodal parameters) 892 CALL CMISSBoundaryConditions_Initialise(BoundaryConditions,Err) 893 CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsSolid,BoundaryConditions,Err) 894 895 !Grab the list of nodes on inner, outer and top surfaces 896 CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_TOP_SURFACE,TopSurfaceNodes,TopNormalXi,Err) 897 CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_INNER_SURFACE,InnerSurfaceNodes,InnerNormalXi,Err) 898 CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_OUTER_SURFACE,OuterSurfaceNodes,OuterNormalXi,Err) 899 900 write(*,*)'TopSurfaceNodes = ',TopSurfaceNodes 901 902 ! ASSIGN BOUNDARY CONDITIONS 903 !Fix base of the ellipsoid in z direction 904 DO NN=1,SIZE(TopSurfaceNodes,1) 905 NODE=TopSurfaceNodes(NN) 906 CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err) 907 IF(NodeDomain==ComputationalNodeNumber) THEN 908 CALL CMISSField_ParameterSetGetNode(GeometricFieldSolid,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, & 909 & ZCoord,Err) 910 CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, & 911 & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err) 912 ENDIF 913 ENDDO 914 915 !Apply inner surface pressure 916 !NOTE: Surface pressure goes into pressure_values_set_type of the DELUDELN type 917 DO NN=1,SIZE(InnerSurfaceNodes,1) 918 NODE=InnerSurfaceNodes(NN) 919 CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err) 920 IF(NodeDomain==ComputationalNodeNumber) THEN 921 CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, & 922 & ABS(InnerNormalXi),CMISS_BOUNDARY_CONDITION_PRESSURE,INNER_PRESSURE,Err) 923 ENDIF 924 ENDDO 925 926 !Apply outer surface pressure 927 DO NN=1,SIZE(OuterSurfaceNodes,1) 928 NODE=OuterSurfaceNodes(NN) 929 CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err) 930 IF(NodeDomain==ComputationalNodeNumber) THEN 931 CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, & 932 & ABS(OuterNor…
Large files files are truncated, but you can click here to view the full file