/FluidMechanics/Stokes/RoutineCheck/Static/src/StaticExample.f90

http://github.com/xyan075/examples · Fortran Modern · 767 lines · 443 code · 102 blank · 222 comment · 0 complexity · 4f55797a4f7249473fdeb08bc16024c9 MD5 · raw file

  1. !> \file
  2. !> \author Sebastian Krittian
  3. !> \brief This is an example program to solve a static Stokes equation using OpenCMISS calls.
  4. !>
  5. !> \section LICENSE
  6. !>
  7. !> Version: MPL 1.1/GPL 2.0/LGPL 2.1
  8. !>
  9. !> The contents of this file are subject to the Mozilla Public License
  10. !> Version 1.1 (the "License"); you may not use this file except in
  11. !> compliance with the License. You may obtain a copy of the License at
  12. !> http://www.mozilla.org/MPL/
  13. !>
  14. !> Software distributed under the License is distributed on an "AS IS"
  15. !> basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
  16. !> License for the specific language governing rights and limitations
  17. !> under the License.
  18. !>
  19. !> The Original Code is OpenCMISS
  20. !>
  21. !> The Initial Developer of the Original Code is University of Auckland,
  22. !> Auckland, New Zealand and University of Oxford, Oxford, United
  23. !> Kingdom. Portions created by the University of Auckland and University
  24. !> of Oxford are Copyright (C) 2007 by the University of Auckland and
  25. !> the University of Oxford. All Rights Reserved.
  26. !>
  27. !> Contributor(s):
  28. !>
  29. !> Alternatively, the contents of this file may be used under the terms of
  30. !> either the GNU General Public License Version 2 or later (the "GPL"), or
  31. !> the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
  32. !> in which case the provisions of the GPL or the LGPL are applicable instead
  33. !> of those above. If you wish to allow use of your version of this file only
  34. !> under the terms of either the GPL or the LGPL, and not to allow others to
  35. !> use your version of this file under the terms of the MPL, indicate your
  36. !> decision by deleting the provisions above and replace them with the notice
  37. !> and other provisions required by the GPL or the LGPL. If you do not delete
  38. !> the provisions above, a recipient may use your version of this file under
  39. !> the terms of any one of the MPL, the GPL or the LGPL.
  40. !>
  41. !> \example FluidMechanics/Stokes/RoutineCheck/Static/src/StaticExample.f90
  42. !! Example program to solve a static Stokes equation using OpenCMISS calls.
  43. !! \htmlinclude FluidMechanics/Stokes/RoutineCheck/Static/history.html
  44. !!
  45. !<
  46. !> Main program
  47. PROGRAM STOKESSTATICEXAMPLE
  48. !
  49. !================================================================================================================================
  50. !
  51. !PROGRAM LIBRARIES
  52. USE OPENCMISS
  53. USE FLUID_MECHANICS_IO_ROUTINES
  54. USE MPI
  55. #ifdef WIN32
  56. USE IFQWINCMISS
  57. #endif
  58. !
  59. !================================================================================================================================
  60. !
  61. !PROGRAM VARIABLES AND TYPES
  62. IMPLICIT NONE
  63. INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumber=1337
  64. TYPE(CMISSFieldType) :: EquationsSetField
  65. !Test program parameters
  66. INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
  67. INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=2
  68. INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=3
  69. INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=4
  70. INTEGER(CMISSIntg), PARAMETER :: GeometricFieldUserNumber=5
  71. INTEGER(CMISSIntg), PARAMETER :: DependentFieldUserNumberStokes=6
  72. INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberStokes=7
  73. INTEGER(CMISSIntg), PARAMETER :: IndependentFieldUserNumberStokes=8
  74. INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberStokes=9
  75. INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=10
  76. INTEGER(CMISSIntg), PARAMETER :: DomainUserNumber=2
  77. INTEGER(CMISSIntg), PARAMETER :: SolverStokesUserNumber=1
  78. INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberStokesMu=1
  79. INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberStokesRho=2
  80. !Program types
  81. TYPE(EXPORT_CONTAINER):: CM
  82. !Program variables
  83. INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
  84. INTEGER(CMISSIntg) :: BASIS_TYPE
  85. INTEGER(CMISSIntg) :: BASIS_NUMBER_SPACE
  86. INTEGER(CMISSIntg) :: BASIS_NUMBER_VELOCITY
  87. INTEGER(CMISSIntg) :: BASIS_NUMBER_PRESSURE
  88. INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_SPACE
  89. INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_VELOCITY
  90. INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_PRESSURE
  91. INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SPACE
  92. INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_VELOCITY
  93. INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_PRESSURE
  94. INTEGER(CMISSIntg) :: MESH_NUMBER_OF_COMPONENTS
  95. INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_SPACE
  96. INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_VELOCITY
  97. INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_PRESSURE
  98. INTEGER(CMISSIntg) :: NUMBER_OF_NODES_SPACE
  99. INTEGER(CMISSIntg) :: NUMBER_OF_NODES_VELOCITY
  100. INTEGER(CMISSIntg) :: NUMBER_OF_NODES_PRESSURE
  101. INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_SPACE
  102. INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_VELOCITY
  103. INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_PRESSURE
  104. INTEGER(CMISSIntg) :: TOTAL_NUMBER_OF_NODES
  105. INTEGER(CMISSIntg) :: TOTAL_NUMBER_OF_ELEMENTS
  106. INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
  107. INTEGER(CMISSIntg) :: RESTART_VALUE
  108. ! INTEGER(CMISSIntg) :: MPI_IERROR
  109. INTEGER(CMISSIntg) :: NUMBER_OF_FIXED_WALL_NODES_STOKES
  110. INTEGER(CMISSIntg) :: NUMBER_OF_INLET_WALL_NODES_STOKES
  111. INTEGER(CMISSIntg) :: EQUATIONS_STOKES_OUTPUT
  112. INTEGER(CMISSIntg) :: COMPONENT_NUMBER
  113. INTEGER(CMISSIntg) :: NODE_NUMBER
  114. INTEGER(CMISSIntg) :: ELEMENT_NUMBER
  115. INTEGER(CMISSIntg) :: NODE_COUNTER
  116. INTEGER(CMISSIntg) :: CONDITION
  117. INTEGER(CMISSIntg) :: LINEAR_SOLVER_STOKES_OUTPUT_TYPE
  118. INTEGER, ALLOCATABLE, DIMENSION(:):: FIXED_WALL_NODES_STOKES
  119. INTEGER, ALLOCATABLE, DIMENSION(:):: INLET_WALL_NODES_STOKES
  120. REAL(CMISSDP) :: INITIAL_FIELD_STOKES(3)
  121. REAL(CMISSDP) :: BOUNDARY_CONDITIONS_STOKES(3)
  122. REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
  123. REAL(CMISSDP) :: RELATIVE_TOLERANCE
  124. REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
  125. REAL(CMISSDP) :: LINESEARCH_ALPHA
  126. REAL(CMISSDP) :: VALUE
  127. REAL(CMISSDP) :: MU_PARAM_STOKES
  128. REAL(CMISSDP) :: RHO_PARAM_STOKES
  129. LOGICAL :: EXPORT_FIELD_IO
  130. LOGICAL :: LINEAR_SOLVER_STOKES_DIRECT_FLAG
  131. LOGICAL :: FIXED_WALL_NODES_STOKES_FLAG
  132. LOGICAL :: INLET_WALL_NODES_STOKES_FLAG
  133. !CMISS variables
  134. !Regions
  135. TYPE(CMISSRegionType) :: Region
  136. TYPE(CMISSRegionType) :: WorldRegion
  137. !Coordinate systems
  138. TYPE(CMISSCoordinateSystemType) :: CoordinateSystem
  139. TYPE(CMISSCoordinateSystemType) :: WorldCoordinateSystem
  140. !Basis
  141. TYPE(CMISSBasisType) :: BasisSpace
  142. TYPE(CMISSBasisType) :: BasisVelocity
  143. TYPE(CMISSBasisType) :: BasisPressure
  144. !Nodes
  145. TYPE(CMISSNodesType) :: Nodes
  146. !Elements
  147. TYPE(CMISSMeshElementsType) :: MeshElementsSpace
  148. TYPE(CMISSMeshElementsType) :: MeshElementsVelocity
  149. TYPE(CMISSMeshElementsType) :: MeshElementsPressure
  150. !Meshes
  151. TYPE(CMISSMeshType) :: Mesh
  152. !Decompositions
  153. TYPE(CMISSDecompositionType) :: Decomposition
  154. !Fields
  155. TYPE(CMISSFieldsType) :: Fields
  156. !Field types
  157. TYPE(CMISSFieldType) :: GeometricField
  158. TYPE(CMISSFieldType) :: DependentFieldStokes
  159. TYPE(CMISSFieldType) :: MaterialsFieldStokes
  160. !Boundary conditions
  161. TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsStokes
  162. !Equations sets
  163. TYPE(CMISSEquationsSetType) :: EquationsSetStokes
  164. !Equations
  165. TYPE(CMISSEquationsType) :: EquationsStokes
  166. !Problems
  167. TYPE(CMISSProblemType) :: Problem
  168. !Control loops
  169. TYPE(CMISSControlLoopType) :: ControlLoop
  170. !Solvers
  171. TYPE(CMISSSolverType) :: LinearSolverStokes
  172. !Solver equations
  173. TYPE(CMISSSolverEquationsType) :: SolverEquationsStokes
  174. #ifdef WIN32
  175. !Quickwin type
  176. LOGICAL :: QUICKWIN_STATUS=.FALSE.
  177. TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
  178. #endif
  179. !Generic CMISS variables
  180. INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber,BoundaryNodeDomain
  181. INTEGER(CMISSIntg) :: EquationsSetIndex
  182. INTEGER(CMISSIntg) :: Err
  183. #ifdef WIN32
  184. !Initialise QuickWin
  185. QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
  186. QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
  187. QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
  188. !Set the window parameters
  189. QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
  190. !If attempt fails set with system estimated values
  191. IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
  192. #endif
  193. !
  194. !================================================================================================================================
  195. !
  196. !INITIALISE OPENCMISS
  197. CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
  198. CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
  199. !
  200. !================================================================================================================================
  201. !
  202. !CHECK COMPUTATIONAL NODE
  203. !Get the computational nodes information
  204. CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
  205. CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
  206. !
  207. !================================================================================================================================
  208. !
  209. !PROBLEM CONTROL PANEL
  210. !Import cmHeart mesh information
  211. CALL FLUID_MECHANICS_IO_READ_CMHEART(CM,Err)
  212. BASIS_NUMBER_SPACE=CM%ID_M
  213. BASIS_NUMBER_VELOCITY=CM%ID_V
  214. BASIS_NUMBER_PRESSURE=CM%ID_P
  215. NUMBER_OF_DIMENSIONS=CM%D
  216. BASIS_TYPE=CM%IT_T
  217. BASIS_XI_INTERPOLATION_SPACE=CM%IT_M
  218. BASIS_XI_INTERPOLATION_VELOCITY=CM%IT_V
  219. BASIS_XI_INTERPOLATION_PRESSURE=CM%IT_P
  220. NUMBER_OF_NODES_SPACE=CM%N_M
  221. NUMBER_OF_NODES_VELOCITY=CM%N_V
  222. NUMBER_OF_NODES_PRESSURE=CM%N_P
  223. TOTAL_NUMBER_OF_NODES=CM%N_T
  224. TOTAL_NUMBER_OF_ELEMENTS=CM%E_T
  225. NUMBER_OF_ELEMENT_NODES_SPACE=CM%EN_M
  226. NUMBER_OF_ELEMENT_NODES_VELOCITY=CM%EN_V
  227. NUMBER_OF_ELEMENT_NODES_PRESSURE=CM%EN_P
  228. !Set initial values
  229. INITIAL_FIELD_STOKES(1)=0.0_CMISSDP
  230. INITIAL_FIELD_STOKES(2)=0.0_CMISSDP
  231. INITIAL_FIELD_STOKES(3)=0.0_CMISSDP
  232. !Set boundary conditions
  233. FIXED_WALL_NODES_STOKES_FLAG=.TRUE.
  234. INLET_WALL_NODES_STOKES_FLAG=.TRUE.
  235. IF(FIXED_WALL_NODES_STOKES_FLAG) THEN
  236. NUMBER_OF_FIXED_WALL_NODES_STOKES=80
  237. ALLOCATE(FIXED_WALL_NODES_STOKES(NUMBER_OF_FIXED_WALL_NODES_STOKES))
  238. FIXED_WALL_NODES_STOKES=(/1,2,3,4,5,7,9,10,11,12,13,14,17,20,24,28,29,30,31,32,33,34,35,37,39, &
  239. & 41,44,46,47,48,50,51,52,53,54,57,60,64,65,66,67,68,70,72,74,76,77,78,79,80,83,86, &
  240. & 89,90,91,92,93,94,95,97,99,101,102,103,104,105,106,107,108,111,114,115,116,117,118, &
  241. & 120,122,123,124,125/)
  242. ENDIF
  243. IF(INLET_WALL_NODES_STOKES_FLAG) THEN
  244. NUMBER_OF_INLET_WALL_NODES_STOKES=9
  245. ALLOCATE(INLET_WALL_NODES_STOKES(NUMBER_OF_INLET_WALL_NODES_STOKES))
  246. INLET_WALL_NODES_STOKES=(/6,15,16,23,36,42,81,82,96/)
  247. !Set initial boundary conditions
  248. BOUNDARY_CONDITIONS_STOKES(1)=0.0_CMISSDP
  249. BOUNDARY_CONDITIONS_STOKES(2)=1.0_CMISSDP
  250. BOUNDARY_CONDITIONS_STOKES(3)=0.0_CMISSDP
  251. ENDIF
  252. !Set material parameters
  253. MU_PARAM_STOKES=1.0_CMISSDP
  254. RHO_PARAM_STOKES=1.0_CMISSDP
  255. !Set interpolation parameters
  256. BASIS_XI_GAUSS_SPACE=3
  257. BASIS_XI_GAUSS_VELOCITY=3
  258. BASIS_XI_GAUSS_PRESSURE=3
  259. !Set output parameter
  260. !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
  261. LINEAR_SOLVER_STOKES_OUTPUT_TYPE=CMISS_SOLVER_NO_OUTPUT
  262. !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
  263. EQUATIONS_STOKES_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
  264. !Set solver parameters
  265. LINEAR_SOLVER_STOKES_DIRECT_FLAG=.FALSE.
  266. RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
  267. ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
  268. DIVERGENCE_TOLERANCE=1.0E20 !default: 1.0E5
  269. MAXIMUM_ITERATIONS=100000 !default: 100000
  270. RESTART_VALUE=3000 !default: 30
  271. LINESEARCH_ALPHA=1.0
  272. !
  273. !================================================================================================================================
  274. !
  275. !COORDINATE SYSTEM
  276. !Start the creation of a new RC coordinate system
  277. CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
  278. CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
  279. !Set the coordinate system dimension
  280. CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NUMBER_OF_DIMENSIONS,Err)
  281. !Finish the creation of the coordinate system
  282. CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
  283. !
  284. !================================================================================================================================
  285. !
  286. !REGION
  287. !Start the creation of a new region
  288. CALL CMISSRegion_Initialise(Region,Err)
  289. CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
  290. !Set the regions coordinate system as defined above
  291. CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
  292. !Finish the creation of the region
  293. CALL CMISSRegion_CreateFinish(Region,Err)
  294. !
  295. !================================================================================================================================
  296. !
  297. !BASES
  298. !Start the creation of new bases
  299. MESH_NUMBER_OF_COMPONENTS=1
  300. CALL CMISSBasis_Initialise(BasisSpace,Err)
  301. CALL CMISSBasis_CreateStart(BASIS_NUMBER_SPACE,BasisSpace,Err)
  302. !Set the basis type (Lagrange/Simplex)
  303. CALL CMISSBasis_TypeSet(BasisSpace,BASIS_TYPE,Err)
  304. !Set the basis xi number
  305. CALL CMISSBasis_NumberOfXiSet(BasisSpace,NUMBER_OF_DIMENSIONS,Err)
  306. !Set the basis xi interpolation and number of Gauss points
  307. IF(NUMBER_OF_DIMENSIONS==2) THEN
  308. CALL CMISSBasis_InterpolationXiSet(BasisSpace,(/BASIS_XI_INTERPOLATION_SPACE,BASIS_XI_INTERPOLATION_SPACE/),Err)
  309. CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisSpace,(/BASIS_XI_GAUSS_SPACE,BASIS_XI_GAUSS_SPACE/),Err)
  310. ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
  311. CALL CMISSBasis_InterpolationXiSet(BasisSpace,(/BASIS_XI_INTERPOLATION_SPACE,BASIS_XI_INTERPOLATION_SPACE, &
  312. & BASIS_XI_INTERPOLATION_SPACE/),Err)
  313. CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisSpace,(/BASIS_XI_GAUSS_SPACE,BASIS_XI_GAUSS_SPACE,BASIS_XI_GAUSS_SPACE/),Err)
  314. ENDIF
  315. !Finish the creation of the basis
  316. CALL CMISSBasis_CreateFinish(BasisSpace,Err)
  317. !Start the creation of another basis
  318. IF(BASIS_XI_INTERPOLATION_VELOCITY==BASIS_XI_INTERPOLATION_SPACE) THEN
  319. BasisVelocity=BasisSpace
  320. ELSE
  321. MESH_NUMBER_OF_COMPONENTS=MESH_NUMBER_OF_COMPONENTS+1
  322. !Initialise a new velocity basis
  323. CALL CMISSBasis_Initialise(BasisVelocity,Err)
  324. !Start the creation of a basis
  325. CALL CMISSBasis_CreateStart(BASIS_NUMBER_VELOCITY,BasisVelocity,Err)
  326. !Set the basis type (Lagrange/Simplex)
  327. CALL CMISSBasis_TypeSet(BasisVelocity,BASIS_TYPE,Err)
  328. !Set the basis xi number
  329. CALL CMISSBasis_NumberOfXiSet(BasisVelocity,NUMBER_OF_DIMENSIONS,Err)
  330. !Set the basis xi interpolation and number of Gauss points
  331. IF(NUMBER_OF_DIMENSIONS==2) THEN
  332. CALL CMISSBasis_InterpolationXiSet(BasisVelocity,(/BASIS_XI_INTERPOLATION_VELOCITY,BASIS_XI_INTERPOLATION_VELOCITY/),Err)
  333. CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisVelocity,(/BASIS_XI_GAUSS_VELOCITY,BASIS_XI_GAUSS_VELOCITY/),Err)
  334. ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
  335. CALL CMISSBasis_InterpolationXiSet(BasisVelocity,(/BASIS_XI_INTERPOLATION_VELOCITY,BASIS_XI_INTERPOLATION_VELOCITY, &
  336. & BASIS_XI_INTERPOLATION_VELOCITY/),Err)
  337. CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisVelocity,(/BASIS_XI_GAUSS_VELOCITY,BASIS_XI_GAUSS_VELOCITY, &
  338. & BASIS_XI_GAUSS_VELOCITY/),Err)
  339. ENDIF
  340. !Finish the creation of the basis
  341. CALL CMISSBasis_CreateFinish(BasisVelocity,Err)
  342. ENDIF
  343. !Start the creation of another basis
  344. IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_SPACE) THEN
  345. BasisPressure=BasisSpace
  346. ELSE IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_VELOCITY) THEN
  347. BasisPressure=BasisVelocity
  348. ELSE
  349. MESH_NUMBER_OF_COMPONENTS=MESH_NUMBER_OF_COMPONENTS+1
  350. !Initialise a new pressure basis
  351. CALL CMISSBasis_Initialise(BasisPressure,Err)
  352. !Start the creation of a basis
  353. CALL CMISSBasis_CreateStart(BASIS_NUMBER_PRESSURE,BasisPressure,Err)
  354. !Set the basis type (Lagrange/Simplex)
  355. CALL CMISSBasis_TypeSet(BasisPressure,BASIS_TYPE,Err)
  356. !Set the basis xi number
  357. CALL CMISSBasis_NumberOfXiSet(BasisPressure,NUMBER_OF_DIMENSIONS,Err)
  358. !Set the basis xi interpolation and number of Gauss points
  359. IF(NUMBER_OF_DIMENSIONS==2) THEN
  360. CALL CMISSBasis_InterpolationXiSet(BasisPressure,(/BASIS_XI_INTERPOLATION_PRESSURE,BASIS_XI_INTERPOLATION_PRESSURE/),Err)
  361. CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisPressure,(/BASIS_XI_GAUSS_PRESSURE,BASIS_XI_GAUSS_PRESSURE/),Err)
  362. ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
  363. CALL CMISSBasis_InterpolationXiSet(BasisPressure,(/BASIS_XI_INTERPOLATION_PRESSURE,BASIS_XI_INTERPOLATION_PRESSURE, &
  364. & BASIS_XI_INTERPOLATION_PRESSURE/),Err)
  365. CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisPressure,(/BASIS_XI_GAUSS_PRESSURE,BASIS_XI_GAUSS_PRESSURE, &
  366. & BASIS_XI_GAUSS_PRESSURE/),Err)
  367. ENDIF
  368. !Finish the creation of the basis
  369. CALL CMISSBasis_CreateFinish(BasisPressure,Err)
  370. ENDIF
  371. !
  372. !================================================================================================================================
  373. !
  374. !MESH
  375. !Start the creation of mesh nodes
  376. CALL CMISSNodes_Initialise(Nodes,Err)
  377. CALL CMISSMesh_Initialise(Mesh,Err)
  378. CALL CMISSNodes_CreateStart(Region,TOTAL_NUMBER_OF_NODES,Nodes,Err)
  379. CALL CMISSNodes_CreateFinish(Nodes,Err)
  380. !Start the creation of the mesh
  381. CALL CMISSMesh_CreateStart(MeshUserNumber,Region,NUMBER_OF_DIMENSIONS,Mesh,Err)
  382. !Set number of mesh elements
  383. CALL CMISSMesh_NumberOfElementsSet(Mesh,TOTAL_NUMBER_OF_ELEMENTS,Err)
  384. !Set number of mesh components
  385. CALL CMISSMesh_NumberOfComponentsSet(Mesh,MESH_NUMBER_OF_COMPONENTS,Err)
  386. !Specify spatial mesh component
  387. CALL CMISSMeshElements_Initialise(MeshElementsSpace,Err)
  388. CALL CMISSMeshElements_Initialise(MeshElementsVelocity,Err)
  389. CALL CMISSMeshElements_Initialise(MeshElementsPressure,Err)
  390. MESH_COMPONENT_NUMBER_SPACE=1
  391. MESH_COMPONENT_NUMBER_VELOCITY=1
  392. MESH_COMPONENT_NUMBER_PRESSURE=1
  393. CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_SPACE,BasisSpace,MeshElementsSpace,Err)
  394. DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
  395. CALL CMISSMeshElements_NodesSet(MeshElementsSpace,ELEMENT_NUMBER,CM%M(ELEMENT_NUMBER,1:NUMBER_OF_ELEMENT_NODES_SPACE),Err)
  396. ENDDO
  397. CALL CMISSMeshElements_CreateFinish(MeshElementsSpace,Err)
  398. !Specify velocity mesh component
  399. IF(BASIS_XI_INTERPOLATION_VELOCITY==BASIS_XI_INTERPOLATION_SPACE) THEN
  400. MeshElementsVelocity=MeshElementsSpace
  401. ELSE
  402. MESH_COMPONENT_NUMBER_VELOCITY=MESH_COMPONENT_NUMBER_SPACE+1
  403. CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_VELOCITY,BasisVelocity,MeshElementsVelocity,Err)
  404. DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
  405. CALL CMISSMeshElements_NodesSet(MeshElementsVelocity,ELEMENT_NUMBER,CM%V(ELEMENT_NUMBER, &
  406. & 1:NUMBER_OF_ELEMENT_NODES_VELOCITY),Err)
  407. ENDDO
  408. CALL CMISSMeshElements_CreateFinish(MeshElementsVelocity,Err)
  409. ENDIF
  410. !Specify pressure mesh component
  411. IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_SPACE) THEN
  412. MeshElementsPressure=MeshElementsSpace
  413. MESH_COMPONENT_NUMBER_PRESSURE=MESH_COMPONENT_NUMBER_SPACE
  414. ELSE IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_VELOCITY) THEN
  415. MeshElementsPressure=MeshElementsVelocity
  416. MESH_COMPONENT_NUMBER_PRESSURE=MESH_COMPONENT_NUMBER_VELOCITY
  417. ELSE
  418. MESH_COMPONENT_NUMBER_PRESSURE=MESH_COMPONENT_NUMBER_VELOCITY+1
  419. CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_PRESSURE,BasisPressure,MeshElementsPressure,Err)
  420. DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
  421. CALL CMISSMeshElements_NodesSet(MeshElementsPressure,ELEMENT_NUMBER,CM%P(ELEMENT_NUMBER, &
  422. & 1:NUMBER_OF_ELEMENT_NODES_PRESSURE),Err)
  423. ENDDO
  424. CALL CMISSMeshElements_CreateFinish(MeshElementsPressure,Err)
  425. ENDIF
  426. !Finish the creation of the mesh
  427. CALL CMISSMesh_CreateFinish(Mesh,Err)
  428. !
  429. !================================================================================================================================
  430. !
  431. !GEOMETRIC FIELD
  432. !Create a decomposition
  433. CALL CMISSDecomposition_Initialise(Decomposition,Err)
  434. CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
  435. !Set the decomposition to be a general decomposition with the specified number of domains
  436. CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
  437. CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfComputationalNodes,Err)
  438. !Finish the decomposition
  439. CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
  440. !Start to create a default (geometric) field on the region
  441. CALL CMISSField_Initialise(GeometricField,Err)
  442. CALL CMISSField_CreateStart(GeometricFieldUserNumber,Region,GeometricField,Err)
  443. !Set the field type
  444. CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
  445. !Set the decomposition to use
  446. CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
  447. !Set the scaling to use
  448. CALL CMISSField_ScalingTypeSet(GeometricField,CMISS_FIELD_NO_SCALING,Err)
  449. !Set the mesh component to be used by the field components.
  450. DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
  451. CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
  452. & MESH_COMPONENT_NUMBER_SPACE,Err)
  453. ENDDO
  454. !Finish creating the field
  455. CALL CMISSField_CreateFinish(GeometricField,Err)
  456. !Update the geometric field parameters
  457. DO NODE_NUMBER=1,NUMBER_OF_NODES_SPACE
  458. DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
  459. VALUE=CM%N(NODE_NUMBER,COMPONENT_NUMBER)
  460. CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE_NUMBER,1,BoundaryNodeDomain,Err)
  461. IF(BoundaryNodeDomain==ComputationalNodeNumber) THEN
  462. CALL CMISSField_ParameterSetUpdateNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
  463. & 1,CMISS_NO_GLOBAL_DERIV,NODE_NUMBER,COMPONENT_NUMBER,VALUE,Err)
  464. ENDIF
  465. ENDDO
  466. ENDDO
  467. CALL CMISSField_ParameterSetUpdateStart(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
  468. CALL CMISSField_ParameterSetUpdateFinish(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
  469. !
  470. !================================================================================================================================
  471. !
  472. !EQUATIONS SETS
  473. !Create the equations set for static Stokes
  474. CALL CMISSEquationsSet_Initialise(EquationsSetStokes,Err)
  475. CALL CMISSField_Initialise(EquationsSetField,Err)
  476. CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberStokes,Region,GeometricField,CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
  477. & CMISS_EQUATIONS_SET_STOKES_EQUATION_TYPE,CMISS_EQUATIONS_SET_STATIC_STOKES_SUBTYPE,EquationsSetFieldUserNumber, &
  478. & EquationsSetField, &
  479. & EquationsSetStokes,Err)
  480. !Set the equations set to be a static Stokes problem
  481. !Finish creating the equations set
  482. CALL CMISSEquationsSet_CreateFinish(EquationsSetStokes,Err)
  483. !
  484. !================================================================================================================================
  485. !
  486. !DEPENDENT FIELDS
  487. !Create the equations set dependent field variables for static Stokes
  488. CALL CMISSField_Initialise(DependentFieldStokes,Err)
  489. CALL CMISSEquationsSet_DependentCreateStart(EquationsSetStokes,DependentFieldUserNumberStokes, &
  490. & DependentFieldStokes,Err)
  491. !Set the mesh component to be used by the field components.
  492. DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
  493. CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
  494. & MESH_COMPONENT_NUMBER_VELOCITY,Err)
  495. CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
  496. & MESH_COMPONENT_NUMBER_VELOCITY,Err)
  497. ENDDO
  498. COMPONENT_NUMBER=NUMBER_OF_DIMENSIONS+1
  499. CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, &
  500. & MESH_COMPONENT_NUMBER_PRESSURE,Err)
  501. CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, &
  502. & MESH_COMPONENT_NUMBER_PRESSURE,Err)
  503. !Finish the equations set dependent field variables
  504. CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetStokes,Err)
  505. !Initialise dependent field
  506. DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
  507. CALL CMISSField_ComponentValuesInitialise(DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
  508. & COMPONENT_NUMBER,INITIAL_FIELD_STOKES(COMPONENT_NUMBER),Err)
  509. ENDDO
  510. !
  511. !================================================================================================================================
  512. !
  513. !MATERIALS FIELDS
  514. !Create the equations set materials field variables for static Stokes
  515. CALL CMISSField_Initialise(MaterialsFieldStokes,Err)
  516. CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetStokes,MaterialsFieldUserNumberStokes, &
  517. & MaterialsFieldStokes,Err)
  518. !Finish the equations set materials field variables
  519. CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetStokes,Err)
  520. CALL CMISSField_ComponentValuesInitialise(MaterialsFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
  521. & MaterialsFieldUserNumberStokesMu,MU_PARAM_STOKES,Err)
  522. CALL CMISSField_ComponentValuesInitialise(MaterialsFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
  523. & MaterialsFieldUserNumberStokesRho,RHO_PARAM_STOKES,Err)
  524. !
  525. !================================================================================================================================
  526. !
  527. !EQUATIONS
  528. !Create the equations set equations
  529. CALL CMISSEquations_Initialise(EquationsStokes,Err)
  530. CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetStokes,EquationsStokes,Err)
  531. !Set the equations matrices sparsity type
  532. CALL CMISSEquations_SparsityTypeSet(EquationsStokes,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
  533. !Set the equations set output
  534. CALL CMISSEquations_OutputTypeSet(EquationsStokes,EQUATIONS_STOKES_OUTPUT,Err)
  535. !Finish the equations set equations
  536. CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetStokes,Err)
  537. !
  538. !================================================================================================================================
  539. !
  540. !PROBLEMS
  541. !Start the creation of a problem.
  542. CALL CMISSProblem_Initialise(Problem,Err)
  543. CALL CMISSControlLoop_Initialise(ControlLoop,Err)
  544. CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
  545. !Set the problem to be a static Stokes problem
  546. CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_FLUID_MECHANICS_CLASS,CMISS_PROBLEM_STOKES_EQUATION_TYPE, &
  547. & CMISS_PROBLEM_STATIC_STOKES_SUBTYPE,Err)
  548. !Finish the creation of a problem.
  549. CALL CMISSProblem_CreateFinish(Problem,Err)
  550. !Start the creation of the problem control loop
  551. CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
  552. !Finish creating the problem control loop
  553. CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
  554. !
  555. !================================================================================================================================
  556. !
  557. !SOLVERS
  558. !Start the creation of the problem solvers
  559. CALL CMISSSolver_Initialise(LinearSolverStokes,Err)
  560. CALL CMISSProblem_SolversCreateStart(Problem,Err)
  561. !Get the linear static solver
  562. CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,SolverStokesUserNumber,LinearSolverStokes,Err)
  563. !Set the output type
  564. CALL CMISSSolver_OutputTypeSet(LinearSolverStokes,LINEAR_SOLVER_STOKES_OUTPUT_TYPE,Err)
  565. !Set the solver settings
  566. IF(LINEAR_SOLVER_STOKES_DIRECT_FLAG) THEN
  567. CALL CMISSSolver_LinearTypeSet(LinearSolverStokes,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
  568. CALL CMISSSolver_LibraryTypeSet(LinearSolverStokes,CMISS_SOLVER_MUMPS_LIBRARY,Err)
  569. ELSE
  570. CALL CMISSSolver_LinearTypeSet(LinearSolverStokes,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
  571. CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverStokes,MAXIMUM_ITERATIONS,Err)
  572. CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverStokes,DIVERGENCE_TOLERANCE,Err)
  573. CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverStokes,RELATIVE_TOLERANCE,Err)
  574. CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverStokes,ABSOLUTE_TOLERANCE,Err)
  575. CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverStokes,RESTART_VALUE,Err)
  576. ENDIF
  577. !Finish the creation of the problem solver
  578. CALL CMISSProblem_SolversCreateFinish(Problem,Err)
  579. !
  580. !================================================================================================================================
  581. !
  582. !SOLVER EQUATIONS
  583. !Start the creation of the problem solver equations
  584. CALL CMISSSolver_Initialise(LinearSolverStokes,Err)
  585. CALL CMISSSolverEquations_Initialise(SolverEquationsStokes,Err)
  586. CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
  587. !Get the linear solver equations
  588. CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,SolverStokesUserNumber,LinearSolverStokes,Err)
  589. CALL CMISSSolver_SolverEquationsGet(LinearSolverStokes,SolverEquationsStokes,Err)
  590. !Set the solver equations sparsity
  591. CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsStokes,CMISS_SOLVER_SPARSE_MATRICES,Err)
  592. !Add in the equations set
  593. CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsStokes,EquationsSetStokes,EquationsSetIndex,Err)
  594. !Finish the creation of the problem solver equations
  595. CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
  596. !
  597. !================================================================================================================================
  598. !
  599. !BOUNDARY CONDITIONS
  600. !Start the creation of the equations set boundary conditions for Stokes
  601. CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsStokes,Err)
  602. CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsStokes,BoundaryConditionsStokes,Err)
  603. !Set fixed wall nodes
  604. IF(FIXED_WALL_NODES_STOKES_FLAG) THEN
  605. DO NODE_COUNTER=1,NUMBER_OF_FIXED_WALL_NODES_STOKES
  606. NODE_NUMBER=FIXED_WALL_NODES_STOKES(NODE_COUNTER)
  607. CONDITION=CMISS_BOUNDARY_CONDITION_FIXED_WALL
  608. CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE_NUMBER,1,BoundaryNodeDomain,Err)
  609. IF(BoundaryNodeDomain==ComputationalNodeNumber) THEN
  610. DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
  611. VALUE=0.0_CMISSDP
  612. CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsStokes,DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,1, &
  613. & CMISS_NO_GLOBAL_DERIV, &
  614. & NODE_NUMBER,COMPONENT_NUMBER,CONDITION,VALUE,Err)
  615. ENDDO
  616. ENDIF
  617. ENDDO
  618. ENDIF
  619. !Set velocity boundary conditions
  620. IF(INLET_WALL_NODES_STOKES_FLAG) THEN
  621. DO NODE_COUNTER=1,NUMBER_OF_INLET_WALL_NODES_STOKES
  622. NODE_NUMBER=INLET_WALL_NODES_STOKES(NODE_COUNTER)
  623. CONDITION=CMISS_BOUNDARY_CONDITION_FIXED_INLET
  624. CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE_NUMBER,1,BoundaryNodeDomain,Err)
  625. IF(BoundaryNodeDomain==ComputationalNodeNumber) THEN
  626. DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
  627. VALUE=BOUNDARY_CONDITIONS_STOKES(COMPONENT_NUMBER)
  628. CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsStokes,DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,1, &
  629. & CMISS_NO_GLOBAL_DERIV, &
  630. & NODE_NUMBER,COMPONENT_NUMBER,CONDITION,VALUE,Err)
  631. ENDDO
  632. ENDIF
  633. ENDDO
  634. ENDIF
  635. !Finish the creation of the equations set boundary conditions
  636. CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsStokes,Err)
  637. !
  638. !================================================================================================================================
  639. !
  640. !RUN SOLVERS
  641. !Turn of PETSc error handling
  642. !CALL PETSC_ERRORHANDLING_SET_ON(ERR,ERROR,*999)
  643. !Solve the problem
  644. WRITE(*,'(A)') "Solving problem..."
  645. CALL CMISSProblem_Solve(Problem,Err)
  646. WRITE(*,'(A)') "Problem solved!"
  647. !
  648. !================================================================================================================================
  649. !
  650. !OUTPUT
  651. EXPORT_FIELD_IO=.TRUE.
  652. IF(EXPORT_FIELD_IO) THEN
  653. WRITE(*,'(A)') "Exporting fields..."
  654. CALL CMISSFields_Initialise(Fields,Err)
  655. CALL CMISSFields_Create(Region,Fields,Err)
  656. CALL CMISSFields_NodesExport(Fields,"StaticStokes","FORTRAN",Err)
  657. CALL CMISSFields_ElementsExport(Fields,"StaticStokes","FORTRAN",Err)
  658. CALL CMISSFields_Finalise(Fields,Err)
  659. WRITE(*,'(A)') "Field exported!"
  660. ENDIF
  661. !Finialise CMISS
  662. CALL CMISSFinalise(Err)
  663. WRITE(*,'(A)') "Program successfully completed."
  664. STOP
  665. END PROGRAM STOKESSTATICEXAMPLE