PageRenderTime 69ms CodeModel.GetById 15ms app.highlight 48ms RepoModel.GetById 1ms app.codeStats 0ms

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