OpenCMISS-Iron Internal API Documentation
Stokes_equations_routines.f90
Go to the documentation of this file.
1 
44 
46 MODULE stokes_equations_routines
47 
48  USE analytic_analysis_routines
49  USE base_routines
50  USE basis_routines
51  USE boundary_conditions_routines
52  USE constants
53  USE control_loop_routines
54  USE distributed_matrix_vector
55  USE domain_mappings
56  USE equations_routines
57  USE equations_mapping_routines
58  USE equations_matrices_routines
59  USE equations_set_constants
60  USE field_routines
61  USE field_io_routines
62  USE fluid_mechanics_io_routines
63  USE input_output
64  USE iso_varying_string
65  USE kinds
66  USE matrix_vector
67  USE node_routines
68  USE problem_constants
69  USE strings
70  USE solver_routines
71  USE timer
72  USE types
73 
74 #include "macros.h"
75 
76  IMPLICIT NONE
77 
78  PRIVATE
79 
80  PUBLIC stokes_equationssetspecificationset
81 
82  PUBLIC stokes_equationssetsolutionmethodset
83 
84  PUBLIC stokes_equations_set_setup
85 
86  PUBLIC stokes_boundaryconditionsanalyticcalculate
87 
88  PUBLIC stokes_problemspecificationset
89 
90  PUBLIC stokes_problem_setup
91 
92  PUBLIC stokes_finite_element_calculate
93 
94  PUBLIC stokes_post_solve
95 
96  PUBLIC stokes_pre_solve
97 
98  PUBLIC stokes_equation_analytic_functions
99 
100 CONTAINS
101 
102 !
103 !================================================================================================================================
104 !
105 
107  SUBROUTINE stokes_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)
108 
109  !Argument variables
110  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
111  INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
112  INTEGER(INTG), INTENT(OUT) :: ERR
113  TYPE(varying_string), INTENT(OUT) :: ERROR
114  !Local Variables
115  TYPE(varying_string) :: LOCAL_ERROR
116 
117  enters("Stokes_EquationsSetSolutionMethodSet",err,error,*999)
118 
119  IF(ASSOCIATED(equations_set)) THEN
120  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
121  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
122  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
123  CALL flagerror("Equations set specification must have 3 entries for a Stokes equations set.", &
124  & err,error,*999)
125  END IF
126  SELECT CASE(equations_set%SPECIFICATION(3))
127  CASE(equations_set_static_stokes_subtype, &
128  & equations_set_laplace_stokes_subtype, &
129  & equations_set_transient_stokes_subtype, &
130  & equations_set_ale_stokes_subtype, &
131  & equations_set_pgm_stokes_subtype)
132  SELECT CASE(solution_method)
133  CASE(equations_set_fem_solution_method)
134  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
135  CASE(equations_set_bem_solution_method)
136  CALL flagerror("Not implemented.",err,error,*999)
137  CASE(equations_set_fd_solution_method)
138  CALL flagerror("Not implemented.",err,error,*999)
139  CASE(equations_set_fv_solution_method)
140  CALL flagerror("Not implemented.",err,error,*999)
141  CASE(equations_set_gfem_solution_method)
142  CALL flagerror("Not implemented.",err,error,*999)
143  CASE(equations_set_gfv_solution_method)
144  CALL flagerror("Not implemented.",err,error,*999)
145  CASE DEFAULT
146  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
147  CALL flagerror(local_error,err,error,*999)
148  END SELECT
149  CASE DEFAULT
150  local_error="The third equations set specification of "// &
151  & trim(numbertovstring(equations_set%SPECIFICATION(3),"*",err,error))// &
152  & " is not valid for a Stokes flow equation of a fluid mechanics equations set."
153  CALL flagerror(local_error,err,error,*999)
154  END SELECT
155  ELSE
156  CALL flagerror("Equations set is not associated.",err,error,*999)
157  ENDIF
158 
159  exits("Stokes_EquationsSetSolutionMethodSet")
160  RETURN
161 999 errorsexits("Stokes_EquationsSetSolutionMethodSet",err,error)
162  RETURN 1
163 
164  END SUBROUTINE stokes_equationssetsolutionmethodset
165 
166 !
167 !================================================================================================================================
168 !
169 
171  SUBROUTINE stokes_equationssetspecificationset(equationsSet,specification,err,error,*)
172 
173  !Argument variables
174  TYPE(equations_set_type), POINTER :: equationsSet
175  INTEGER(INTG), INTENT(IN) :: specification(:)
176  INTEGER(INTG), INTENT(OUT) :: err
177  TYPE(varying_string), INTENT(OUT) :: error
178  !Local Variables
179  TYPE(varying_string) :: localError
180  INTEGER(INTG) :: subtype
181 
182  enters("Stokes_EquationsSetSpecificationSet",err,error,*999)
183 
184  IF(ASSOCIATED(equationsset)) THEN
185  IF(SIZE(specification,1)/=3) THEN
186  CALL flagerror("Equations set specification must have three entries for a Stokes type equations set.", &
187  & err,error,*999)
188  END IF
189  subtype=specification(3)
190  SELECT CASE(subtype)
191  CASE(equations_set_static_stokes_subtype, &
192  & equations_set_laplace_stokes_subtype, &
193  & equations_set_transient_stokes_subtype, &
194  & equations_set_ale_stokes_subtype, &
195  & equations_set_pgm_stokes_subtype)
196  !ok
197  CASE(equations_set_optimised_stokes_subtype)
198  CALL flagerror("Not implemented yet.",err,error,*999)
199  CASE DEFAULT
200  localerror="The third equations set specification of "//trim(numbertovstring(specification(3),"*",err,error))// &
201  & " is not valid for Stokes flow of a fluid mechanics equations set."
202  CALL flagerror(localerror,err,error,*999)
203  END SELECT
204  !Set full specification
205  IF(ALLOCATED(equationsset%specification)) THEN
206  CALL flagerror("Equations set specification is already allocated.",err,error,*999)
207  ELSE
208  ALLOCATE(equationsset%specification(3),stat=err)
209  IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
210  END IF
211  equationsset%specification(1:3)=[equations_set_fluid_mechanics_class,equations_set_stokes_equation_type,subtype]
212  ELSE
213  CALL flagerror("Equations set is not associated.",err,error,*999)
214  END IF
215 
216  exits("Stokes_EquationsSetSpecificationSet")
217  RETURN
218 999 errors("Stokes_EquationsSetSpecificationSet",err,error)
219  exits("Stokes_EquationsSetSpecificationSet")
220  RETURN 1
221 
222  END SUBROUTINE stokes_equationssetspecificationset
223 
224 !
225 !================================================================================================================================
226 !
227 
229  SUBROUTINE stokes_equations_set_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
230 
231  !Argument variables
232  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
233  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
234  INTEGER(INTG), INTENT(OUT) :: ERR
235  TYPE(varying_string), INTENT(OUT) :: ERROR
236  !Local Variables
237  INTEGER(INTG) :: GEOMETRIC_SCALING_TYPE,GEOMETRIC_MESH_COMPONENT
238  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
239  TYPE(equations_type), POINTER :: EQUATIONS
240  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
241  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
242  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
243  TYPE(varying_string) :: LOCAL_ERROR
244  INTEGER(INTG):: DEPENDENT_FIELD_NUMBER_OF_VARIABLES,DEPENDENT_FIELD_NUMBER_OF_COMPONENTS
245  INTEGER(INTG):: INDEPENDENT_FIELD_NUMBER_OF_VARIABLES,INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS
246  INTEGER(INTG):: NUMBER_OF_DIMENSIONS,GEOMETRIC_COMPONENT_NUMBER
247  INTEGER(INTG):: MATERIAL_FIELD_NUMBER_OF_VARIABLES,MATERIAL_FIELD_NUMBER_OF_COMPONENTS,I
248 
249  enters("STOKES_EQUATION_SET_SETUP",err,error,*999)
250 
251  NULLIFY(equations)
252  NULLIFY(equations_mapping)
253  NULLIFY(equations_matrices)
254  NULLIFY(geometric_decomposition)
255 
256  IF(ASSOCIATED(equations_set)) THEN
257  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
258  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
259  ELSE IF(SIZE(equations_set%SPECIFICATION,1)<3) THEN
260  CALL flagerror("Equations set specification must have >= 3 entries for a Stokes flow equations set.", &
261  & err,error,*999)
262  END IF
263  SELECT CASE(equations_set%SPECIFICATION(3))
264  !Select Stokes subtypes
265  CASE(equations_set_static_stokes_subtype, &
266  & equations_set_transient_stokes_subtype, &
267  & equations_set_laplace_stokes_subtype, &
268  & equations_set_ale_stokes_subtype, &
269  & equations_set_pgm_stokes_subtype)
270  SELECT CASE(equations_set_setup%SETUP_TYPE)
271  !Set solution method
272  CASE(equations_set_setup_initial_type)
273  SELECT CASE(equations_set%SPECIFICATION(3))
274  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype, &
275  & equations_set_transient_stokes_subtype,equations_set_ale_stokes_subtype, &
276  & equations_set_pgm_stokes_subtype)
277  SELECT CASE(equations_set_setup%ACTION_TYPE)
278  CASE(equations_set_setup_start_action)
279  CALL stokes_equationssetsolutionmethodset(equations_set ,&
280  & equations_set_fem_solution_method,err,error,*999)
281  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
282  CASE(equations_set_setup_finish_action)
283  !!TODO: Check valid setup
284  CASE DEFAULT
285  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE, &
286  & "*",err,error))// " for a setup type of "//trim(number_to_vstring(equations_set_setup% &
287  & setup_type,"*",err,error))// " is invalid for a standard Stokes fluid."
288  CALL flagerror(local_error,err,error,*999)
289  END SELECT
290  CASE DEFAULT
291  local_error="The third equations set specification of "// &
292  & trim(numbertovstring(equations_set%SPECIFICATION(3),"*", &
293  & err,error))//" is invalid for a Stokes flow equations set."
294  CALL flagerror(local_error,err,error,*999)
295  END SELECT
296  CASE(equations_set_setup_geometry_type)
297  !Set geometric field
298  SELECT CASE(equations_set%SPECIFICATION(3))
299  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype, &
300  & equations_set_transient_stokes_subtype,equations_set_ale_stokes_subtype, &
301  & equations_set_pgm_stokes_subtype)
302  !Do nothing???
303  CASE DEFAULT
304  local_error="The third equations set specification of "// &
305  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*", &
306  & err,error))//" is invalid for a Stokes flow equations set."
307  CALL flagerror(local_error,err,error,*999)
308  END SELECT
309  CASE(equations_set_setup_dependent_type)
310  !Set dependent field
311  SELECT CASE(equations_set%SPECIFICATION(3))
312  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype, &
313  & equations_set_transient_stokes_subtype,equations_set_ale_stokes_subtype, &
314  & equations_set_pgm_stokes_subtype)
315  SELECT CASE(equations_set_setup%ACTION_TYPE)
316  !Set start action
317  CASE(equations_set_setup_start_action)
318  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
319  !Create the auto created dependent field
320  !start field creation with name 'DEPENDENT_FIELD'
321  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
322  & equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
323  !start creation of a new field
324  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
325  !label the field
326  CALL field_label_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,"Dependent Field",err,error,*999)
327  !define new created field to be dependent
328  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
329  & field_dependent_type,err,error,*999)
330  !look for decomposition rule already defined
331  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
332  & err,error,*999)
333  !apply decomposition rule found on new created field
334  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
335  & geometric_decomposition,err,error,*999)
336  !point new field to geometric field
337  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
338  & geometric_field,err,error,*999)
339  !set number of variables to 2 (1 for U and one for DELUDELN)
340  dependent_field_number_of_variables=2
341  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
342  & dependent_field_number_of_variables,err,error,*999)
343  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,[field_u_variable_type, &
344  & field_deludeln_variable_type],err,error,*999)
345  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
346  & field_vector_dimension_type,err,error,*999)
347  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
348  & field_vector_dimension_type,err,error,*999)
349  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
350  & field_dp_type,err,error,*999)
351  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
352  & field_dp_type,err,error,*999)
353  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
354  & number_of_dimensions,err,error,*999)
355  !calculate number of components with one component for each dimension and one for pressure
356  dependent_field_number_of_components=number_of_dimensions+1
357  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
358  & field_u_variable_type,dependent_field_number_of_components,err,error,*999)
359  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
360  & field_deludeln_variable_type,dependent_field_number_of_components,err,error,*999)
361  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
362  & 1,geometric_mesh_component,err,error,*999)
363  !Default to the geometric interpolation setup
364  DO i=1,dependent_field_number_of_components
365  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD, &
366  & field_u_variable_type,i,geometric_mesh_component,err,error,*999)
367  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD, &
368  & field_deludeln_variable_type,i,geometric_mesh_component,err,error,*999)
369  END DO
370  SELECT CASE(equations_set%SOLUTION_METHOD)
371  !Specify fem solution method
372  CASE(equations_set_fem_solution_method)
373  DO i=1,dependent_field_number_of_components
374  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
375  & field_u_variable_type,i,field_node_based_interpolation,err,error,*999)
376  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
377  & field_deludeln_variable_type,i,field_node_based_interpolation,err,error,*999)
378  END DO
379  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
380  & err,error,*999)
381  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type, &
382  & err,error,*999)
383  !Other solutions not defined yet
384  CASE DEFAULT
385  local_error="The solution method of " &
386  & //trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// " is invalid."
387  CALL flagerror(local_error,err,error,*999)
388  END SELECT
389  ELSE
390  !Check the user specified field
391  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
392  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
393  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
394  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type, &
395  & field_deludeln_variable_type],err,error,*999)
396  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
397  & err,error,*999)
398  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
399  & field_vector_dimension_type,err,error,*999)
400  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
401  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type, &
402  & err,error,*999)
403  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
404  & number_of_dimensions,err,error,*999)
405  !calculate number of components with one component for each dimension and one for pressure
406  dependent_field_number_of_components=number_of_dimensions+1
407  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
408  & dependent_field_number_of_components,err,error,*999)
409  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
410  & dependent_field_number_of_components,err,error,*999)
411  SELECT CASE(equations_set%SOLUTION_METHOD)
412  CASE(equations_set_fem_solution_method)
413  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
414  & field_node_based_interpolation,err,error,*999)
415  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
416  & field_node_based_interpolation,err,error,*999)
417  CASE DEFAULT
418  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
419  &"*",err,error))//" is invalid."
420  CALL flagerror(local_error,err,error,*999)
421  END SELECT
422  ENDIF
423  !Specify finish action
424  CASE(equations_set_setup_finish_action)
425  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
426  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
427  ENDIF
428  CASE DEFAULT
429  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
430  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
431  & " is invalid for a standard Stokes fluid"
432  CALL flagerror(local_error,err,error,*999)
433  END SELECT
434  CASE DEFAULT
435  local_error="The third equations set specification of "// &
436  & trim(numbertovstring(equations_set%SPECIFICATION(3),"*", &
437  & err,error))//" is invalid for a Stokes flow equations set."
438  CALL flag_error(local_error,err,error,*999)
439  END SELECT
440  CASE(equations_set_setup_independent_type)
441  !define an independent field for ALE information
442  SELECT CASE(equations_set%SPECIFICATION(3))
443  CASE(equations_set_ale_stokes_subtype,equations_set_pgm_stokes_subtype)
444  SELECT CASE(equations_set_setup%ACTION_TYPE)
445  !Set start action
446  CASE(equations_set_setup_start_action)
447  IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
448  !Create the auto created independent field
449  !start field creation with name 'INDEPENDENT_FIELD'
450  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
451  & equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
452  !start creation of a new field
453  CALL field_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_general_type,err,error,*999)
454  !label the field
455  CALL field_label_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,"Independent Field",err,error,*999)
456  !define new created field to be independent
457  CALL field_dependent_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
458  & field_independent_type,err,error,*999)
459  !look for decomposition rule already defined
460  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
461  & err,error,*999)
462  !apply decomposition rule found on new created field
463  CALL field_mesh_decomposition_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
464  & geometric_decomposition,err,error,*999)
465  !point new field to geometric field
466  CALL field_geometric_field_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,equations_set% &
467  & geometry%GEOMETRIC_FIELD,err,error,*999)
468  !set number of variables to 1 (1 for U)
469  independent_field_number_of_variables=1
470  CALL field_number_of_variables_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
471  & independent_field_number_of_variables,err,error,*999)
472  CALL field_variable_types_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
473  & [field_u_variable_type],err,error,*999)
474  CALL field_dimension_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
475  & field_vector_dimension_type,err,error,*999)
476  CALL field_data_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
477  & field_dp_type,err,error,*999)
478  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
479  & number_of_dimensions,err,error,*999)
480  !calculate number of components with one component for each dimension
481  independent_field_number_of_components=number_of_dimensions
482  CALL field_number_of_components_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
483  & field_u_variable_type,independent_field_number_of_components,err,error,*999)
484  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
485  & 1,geometric_mesh_component,err,error,*999)
486  !Default to the geometric interpolation setup
487  DO i=1,independent_field_number_of_components
488  CALL field_component_mesh_component_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
489  & field_u_variable_type,i,geometric_mesh_component,err,error,*999)
490  END DO
491  SELECT CASE(equations_set%SOLUTION_METHOD)
492  !Specify fem solution method
493  CASE(equations_set_fem_solution_method)
494  DO i=1,independent_field_number_of_components
495  CALL field_component_interpolation_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
496  & field_u_variable_type,i,field_node_based_interpolation,err,error,*999)
497  END DO
498  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
499  & err,error,*999)
500  CALL field_scaling_type_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,geometric_scaling_type, &
501  & err,error,*999)
502  !Other solutions not defined yet
503  CASE DEFAULT
504  local_error="The solution method of " &
505  & //trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// " is invalid."
506  CALL flagerror(local_error,err,error,*999)
507  END SELECT
508  ELSE
509  !Check the user specified field
510  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
511  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
512  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
513  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
514  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
515  & err,error,*999)
516  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
517  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
518  & number_of_dimensions,err,error,*999)
519  !calculate number of components with one component for each dimension and one for pressure
520  independent_field_number_of_components=number_of_dimensions
521  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
522  & independent_field_number_of_components,err,error,*999)
523  SELECT CASE(equations_set%SOLUTION_METHOD)
524  CASE(equations_set_fem_solution_method)
525  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
526  & field_node_based_interpolation,err,error,*999)
527  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
528  & field_node_based_interpolation,err,error,*999)
529  CASE DEFAULT
530  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
531  &"*",err,error))//" is invalid."
532  CALL flagerror(local_error,err,error,*999)
533  END SELECT
534  ENDIF
535  !Specify finish action
536  CASE(equations_set_setup_finish_action)
537  IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
538  CALL field_create_finish(equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
539  CALL field_parameter_set_create(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
540  & field_mesh_displacement_set_type,err,error,*999)
541  CALL field_parameter_set_create(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
542  & field_mesh_velocity_set_type,err,error,*999)
543  CALL field_parameter_set_create(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
544  & field_boundary_set_type,err,error,*999)
545  ENDIF
546  CASE DEFAULT
547  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
548  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
549  & " is invalid for a standard Stokes fluid"
550  CALL flagerror(local_error,err,error,*999)
551  END SELECT
552  CASE DEFAULT
553  local_error="The third equations set specification of "// &
554  & trim(numbertovstring(equations_set%SPECIFICATION(3),"*", &
555  & err,error))//" is invalid for a Stokes flow equations set."
556  CALL flagerror(local_error,err,error,*999)
557  END SELECT
558  !Define analytic part for Stokes problem
559  CASE(equations_set_setup_analytic_type)
560  SELECT CASE(equations_set%SPECIFICATION(3))
561  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype, &
562  & equations_set_transient_stokes_subtype)
563  SELECT CASE(equations_set_setup%ACTION_TYPE)
564  !Set start action
565  CASE(equations_set_setup_start_action)
566  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
567  IF(ASSOCIATED(equations_set%DEPENDENT%DEPENDENT_FIELD)) THEN
568  IF(ASSOCIATED(equations_set%GEOMETRY%GEOMETRIC_FIELD)) THEN
569  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
570  & number_of_dimensions,err,error,*999)
571  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
572  CASE(equations_set_stokes_equation_two_dim_1)
573  !Set analtyic function type
574  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_two_dim_1
575  CASE(equations_set_stokes_equation_two_dim_2)
576  !Set analtyic function type
577  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_two_dim_2
578  CASE(equations_set_stokes_equation_two_dim_3)
579  !Set analtyic function type
580  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_two_dim_3
581  CASE(equations_set_stokes_equation_two_dim_4)
582  !Set analtyic function type
583  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_two_dim_4
584  CASE(equations_set_stokes_equation_two_dim_5)
585  !Set analtyic function type
586  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_two_dim_5
587  CASE(equations_set_stokes_equation_three_dim_1)
588  !Set analtyic function type
589  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_three_dim_1
590  CASE(equations_set_stokes_equation_three_dim_2)
591  !Set analtyic function type
592  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_three_dim_2
593  CASE(equations_set_stokes_equation_three_dim_3)
594  !Set analtyic function type
595  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_three_dim_3
596  CASE(equations_set_stokes_equation_three_dim_4)
597  !Set analtyic function type
598  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_three_dim_4
599  CASE(equations_set_stokes_equation_three_dim_5)
600  !Set analtyic function type
601  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_stokes_equation_three_dim_5
602  CASE DEFAULT
603  local_error="The specified analytic function type of "// &
604  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
605  & " is invalid for an analytic Stokes problem."
606  CALL flagerror(local_error,err,error,*999)
607  END SELECT
608  ELSE
609  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
610  ENDIF
611  ELSE
612  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
613  ENDIF
614  ELSE
615  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
616  ENDIF
617  CASE(equations_set_setup_finish_action)
618  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
619  IF(ASSOCIATED(equations_set%ANALYTIC%ANALYTIC_FIELD)) THEN
620  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
621  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
622  ENDIF
623  ENDIF
624  ELSE
625  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
626  ENDIF
627  CASE DEFAULT
628  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
629  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
630  & " is invalid for an analytic Stokes problem."
631  CALL flagerror(local_error,err,error,*999)
632  END SELECT
633  CASE DEFAULT
634  local_error="The third equations set specification of "// &
635  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
636  " is invalid for a Stokes flow equations set."
637  CALL flag_error(local_error,err,error,*999)
638  END SELECT
639  !Define materials field
640  CASE(equations_set_setup_materials_type)
641  SELECT CASE(equations_set%SPECIFICATION(3))
642  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype, &
643  & equations_set_transient_stokes_subtype,equations_set_ale_stokes_subtype, &
644  & equations_set_pgm_stokes_subtype)
645  !variable X with has Y components, here Y represents viscosity only
646  material_field_number_of_variables=1!X
647  material_field_number_of_components=2!Y
648  SELECT CASE(equations_set_setup%ACTION_TYPE)
649  !Specify start action
650  CASE(equations_set_setup_start_action)
651  equations_materials=>equations_set%MATERIALS
652  IF(ASSOCIATED(equations_materials)) THEN
653  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
654  !Create the auto created materials field
655  !start field creation with name 'MATERIAL_FIELD'
656  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set% &
657  & materials%MATERIALS_FIELD,err,error,*999)
658  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
659  !label the field
660  CALL field_label_set_and_lock(equations_materials%MATERIALS_FIELD,"Materials Field",err,error,*999)
661  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type, &
662  & err,error,*999)
663  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
664  & err,error,*999)
665  !apply decomposition rule found on new created field
666  CALL field_mesh_decomposition_set_and_lock(equations_set%MATERIALS%MATERIALS_FIELD, &
667  & geometric_decomposition,err,error,*999)
668  !point new field to geometric field
669  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
670  & geometric_field,err,error,*999)
671  CALL field_number_of_variables_set(equations_materials%MATERIALS_FIELD, &
672  & material_field_number_of_variables,err,error,*999)
673  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type], &
674  & err,error,*999)
675  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
676  & field_vector_dimension_type,err,error,*999)
677  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
678  & field_dp_type,err,error,*999)
679  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
680  & material_field_number_of_components,err,error,*999)
681  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
682  & 1,geometric_component_number,err,error,*999)
683  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
684  & 1,geometric_component_number,err,error,*999)
685  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
686  & 1,field_constant_interpolation,err,error,*999)
687  !Default the field scaling to that of the geometric field
688  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
689  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
690  ELSE
691  !Check the user specified field
692  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
693  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
694  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
695  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
696  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
697  & err,error,*999)
698  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
699  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
700  & number_of_dimensions,err,error,*999)
701  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
702  ENDIF
703  ELSE
704  CALL flagerror("Equations set materials is not associated.",err,error,*999)
705  END IF
706  !Specify start action
707  CASE(equations_set_setup_finish_action)
708  equations_materials=>equations_set%MATERIALS
709  IF(ASSOCIATED(equations_materials)) THEN
710  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
711  !Finish creating the materials field
712  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
713  !Set the default values for the materials field
714  !First set the mu values to 0.001
715  !MATERIAL_FIELD_NUMBER_OF_COMPONENTS
716  ! viscosity=1
717  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
718  & field_values_set_type,1,1.0_dp,err,error,*999)
719  ! density=2
720 !\todo: Initialise fields properly
721  ! CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, &
722  ! & FIELD_VALUES_SET_TYPE,2,100.0_DP,ERR,ERROR,*999)
723  ENDIF
724  ELSE
725  CALL flagerror("Equations set materials is not associated.",err,error,*999)
726  ENDIF
727  CASE DEFAULT
728  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
729  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
730  & " is invalid for Stokes equation."
731  CALL flagerror(local_error,err,error,*999)
732  END SELECT
733  CASE DEFAULT
734  local_error="The third equations set specification of "// &
735  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
736  & " is invalid for a Stokes flow equations set."
737  CALL flagerror(local_error,err,error,*999)
738  END SELECT
739  !Define the source field
740  CASE(equations_set_setup_source_type)
741  SELECT CASE(equations_set%SPECIFICATION(3))
742  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype, &
743  & equations_set_transient_stokes_subtype,equations_set_ale_stokes_subtype, &
744  & equations_set_pgm_stokes_subtype)
745  !TO DO: INCLUDE GRAVITY AS SOURCE TYPE
746  SELECT CASE(equations_set_setup%ACTION_TYPE)
747  CASE(equations_set_setup_start_action)
748  !Do nothing
749  CASE(equations_set_setup_finish_action)
750  !Do nothing
751  !? Maybe set finished flag????
752  CASE DEFAULT
753  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
754  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
755  & " is invalid for a standard Stokes fluid."
756  CALL flagerror(local_error,err,error,*999)
757  END SELECT
758  CASE DEFAULT
759  local_error="The third equations set specification of "// &
760  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
761  & " is invalid for a Stokes flow equations set."
762  CALL flag_error(local_error,err,error,*999)
763  END SELECT
764  !Define equations type
765  CASE(equations_set_setup_equations_type)
766  SELECT CASE(equations_set%SPECIFICATION(3))
767  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype)
768  SELECT CASE(equations_set_setup%ACTION_TYPE)
769  CASE(equations_set_setup_start_action)
770  equations_materials=>equations_set%MATERIALS
771  IF(ASSOCIATED(equations_materials)) THEN
772  IF(equations_materials%MATERIALS_FINISHED) THEN
773  CALL equations_create_start(equations_set,equations,err,error,*999)
774  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
775  CALL equations_time_dependence_type_set(equations,equations_static,err,error,*999)
776  ELSE
777  CALL flagerror("Equations set materials has not been finished.",err,error,*999)
778  ENDIF
779  ELSE
780  CALL flagerror("Equations materials is not associated.",err,error,*999)
781  ENDIF
782  CASE(equations_set_setup_finish_action)
783  SELECT CASE(equations_set%SOLUTION_METHOD)
784  CASE(equations_set_fem_solution_method)
785  !Finish the equations creation
786  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
787  CALL equations_create_finish(equations,err,error,*999)
788  !Create the equations mapping.
789  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
790  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
791  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,[field_u_variable_type], &
792  & err,error,*999)
793  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type, &
794  & err,error,*999)
795  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
796  !Create the equations matrices
797  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
798  SELECT CASE(equations%SPARSITY_TYPE)
799  CASE(equations_matrices_full_matrices)
800  CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_block_storage_type], &
801  & err,error,*999)
802  CASE(equations_matrices_sparse_matrices)
803  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
804  & [matrix_compressed_row_storage_type],err,error,*999)
805  CALL equationsmatrices_linearstructuretypeset(equations_matrices, &
806  & [equations_matrix_fem_structure],err,error,*999)
807  CASE DEFAULT
808  local_error="The equations matrices sparsity type of "// &
809  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
810  CALL flagerror(local_error,err,error,*999)
811  END SELECT
812  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
813  CASE DEFAULT
814  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
815  & "*",err,error))//" is invalid."
816  CALL flagerror(local_error,err,error,*999)
817  END SELECT
818  CASE DEFAULT
819  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
820  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
821  & " is invalid for a steady Laplace equation."
822  CALL flagerror(local_error,err,error,*999)
823  END SELECT
824  CASE(equations_set_transient_stokes_subtype,equations_set_ale_stokes_subtype,equations_set_pgm_stokes_subtype)
825  SELECT CASE(equations_set_setup%ACTION_TYPE)
826  CASE(equations_set_setup_start_action)
827  equations_materials=>equations_set%MATERIALS
828  IF(ASSOCIATED(equations_materials)) THEN
829  IF(equations_materials%MATERIALS_FINISHED) THEN
830  CALL equations_create_start(equations_set,equations,err,error,*999)
831  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
832  CALL equations_time_dependence_type_set(equations,equations_first_order_dynamic,err,error,*999)
833  ELSE
834  CALL flagerror("Equations set materials has not been finished.",err,error,*999)
835  ENDIF
836  ELSE
837  CALL flagerror("Equations materials is not associated.",err,error,*999)
838  ENDIF
839  CASE(equations_set_setup_finish_action)
840  SELECT CASE(equations_set%SOLUTION_METHOD)
841  CASE(equations_set_fem_solution_method)
842  !Finish the equations creation
843  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
844  CALL equations_create_finish(equations,err,error,*999)
845  !Create the equations mapping.
846  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
847  CALL equations_mapping_dynamic_matrices_set(equations_mapping,.true.,.true.,err,error,*999)
848  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
849  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type, &
850  & err,error,*999)
851  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
852  !Create the equations matrices
853  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
854  !Set up matrix storage and structure
855  IF(equations%LUMPING_TYPE==equations_lumped_matrices) THEN
856  !Set up lumping
857  CALL equations_matrices_dynamic_lumping_type_set(equations_matrices, &
858  & [equations_matrix_unlumped,equations_matrix_lumped],err,error,*999)
859  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
860  & [distributed_matrix_compressed_row_storage_type,distributed_matrix_diagonal_storage_type] &
861  & ,err,error,*999)
862  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
863  & [equations_matrix_fem_structure,equations_matrix_diagonal_structure],err,error,*999)
864  ELSE
865  SELECT CASE(equations%SPARSITY_TYPE)
866  CASE(equations_matrices_full_matrices)
867  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
868  & [distributed_matrix_block_storage_type,distributed_matrix_block_storage_type],err,error,*999)
869  CASE(equations_matrices_sparse_matrices)
870  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
871  & [distributed_matrix_compressed_row_storage_type, &
872  & distributed_matrix_compressed_row_storage_type],err,error,*999)
873  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
874  & [equations_matrix_fem_structure,equations_matrix_fem_structure],err,error,*999)
875  CASE DEFAULT
876  local_error="The equations matrices sparsity type of "// &
877  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
878  CALL flagerror(local_error,err,error,*999)
879  END SELECT
880  ENDIF
881  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
882  CASE DEFAULT
883  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*", &
884  & err,error))//" is invalid."
885  CALL flagerror(local_error,err,error,*999)
886  END SELECT
887  CASE DEFAULT
888  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
889  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
890  & " is invalid for a Stokes equation."
891  CALL flagerror(local_error,err,error,*999)
892  END SELECT
893  CASE DEFAULT
894  local_error="The third equations set specification of "// &
895  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
896  & " is invalid for a Stokes flow equations set."
897  CALL flagerror(local_error,err,error,*999)
898  END SELECT
899  CASE DEFAULT
900  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
901  & " is invalid for a standard Stokes fluid."
902  CALL flagerror(local_error,err,error,*999)
903  END SELECT
904  CASE DEFAULT
905  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
906  & " does not equal a standard Stokes fluid subtype."
907  CALL flagerror(local_error,err,error,*999)
908  END SELECT
909  ELSE
910  CALL flagerror("Equations set is not associated.",err,error,*999)
911  ENDIF
912 
913  exits("STOKES_EQUATIONS_SET_SETUP")
914  RETURN
915 999 errorsexits("STOKES_EQUATIONS_SET_SETUP",err,error)
916  RETURN 1
917  END SUBROUTINE stokes_equations_set_setup
918 
919 !
920 !================================================================================================================================
921 !
922 
924  SUBROUTINE stokes_problemspecificationset(problem,problemSpecification,err,error,*)
925 
926  !Argument variables
927  TYPE(problem_type), POINTER :: problem
928  INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
929  INTEGER(INTG), INTENT(OUT) :: err
930  TYPE(varying_string), INTENT(OUT) :: error
931  !Local Variables
932  TYPE(varying_string) :: localError
933  INTEGER(INTG) :: problemSubtype
934 
935  enters("Stokes_ProblemSpecificationSet",err,error,*999)
936 
937  IF(ASSOCIATED(problem)) THEN
938  IF(SIZE(problemspecification,1)==3) THEN
939  problemsubtype=problemspecification(3)
940  SELECT CASE(problemsubtype)
941  CASE(problem_static_stokes_subtype, &
942  & problem_laplace_stokes_subtype, &
943  & problem_transient_stokes_subtype, &
944  & problem_ale_stokes_subtype, &
945  & problem_pgm_stokes_subtype)
946  !All ok
947  CASE(problem_optimised_stokes_subtype)
948  CALL flagerror("Not implemented yet.",err,error,*999)
949  CASE DEFAULT
950  localerror="The third problem specification of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
951  & " is not valid for a Stokes flow fluid mechanics problem."
952  CALL flagerror(localerror,err,error,*999)
953  END SELECT
954  IF(ALLOCATED(problem%specification)) THEN
955  CALL flagerror("Problem specification is already allocated.",err,error,*999)
956  ELSE
957  ALLOCATE(problem%specification(3),stat=err)
958  IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
959  ENDIF
960  problem%specification(1:3)=[problem_fluid_mechanics_class,problem_stokes_equation_type,problemsubtype]
961  ELSE
962  CALL flagerror("Stokes flow problem specification must have >= 3 entries.",err,error,*999)
963  ENDIF
964  ELSE
965  CALL flagerror("Problem is not associated.",err,error,*999)
966  ENDIF
967 
968  exits("Stokes_ProblemSpecificationSet")
969  RETURN
970 999 errorsexits("Stokes_ProblemSpecificationSet",err,error)
971  RETURN 1
972 
973  END SUBROUTINE stokes_problemspecificationset
974 
975 !
976 !================================================================================================================================
977 !
978 
980  SUBROUTINE stokes_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
981 
982  !Argument variables
983  TYPE(problem_type), POINTER :: PROBLEM
984  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
985  INTEGER(INTG), INTENT(OUT) :: ERR
986  TYPE(varying_string), INTENT(OUT) :: ERROR
987  !Local Variables
988  TYPE(varying_string) :: LOCAL_ERROR
989  TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
990  TYPE(solver_type), POINTER :: SOLVER, MESH_SOLVER
991  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS,MESH_SOLVER_EQUATIONS
992  TYPE(solvers_type), POINTER :: SOLVERS
993 
994  enters("STOKES_PROBLEM_SETUP",err,error,*999)
995 
996  NULLIFY(control_loop)
997  NULLIFY(solver)
998  NULLIFY(mesh_solver)
999  NULLIFY(solver_equations)
1000  NULLIFY(mesh_solver_equations)
1001  NULLIFY(solvers)
1002  IF(ASSOCIATED(problem)) THEN
1003  IF(.NOT.ALLOCATED(problem%specification)) THEN
1004  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1005  ELSE IF(SIZE(problem%specification,1)<3) THEN
1006  CALL flagerror("Problem specification must have three entries for a Stokes problem.",err,error,*999)
1007  END IF
1008  SELECT CASE(problem%SPECIFICATION(3))
1009  !Set problem subtype for steady state Stokes problems
1010  CASE(problem_static_stokes_subtype,problem_laplace_stokes_subtype)
1011  SELECT CASE(problem_setup%SETUP_TYPE)
1012  CASE(problem_setup_initial_type)
1013  SELECT CASE(problem_setup%ACTION_TYPE)
1014  CASE(problem_setup_start_action)
1015  !Do nothing????
1016  CASE(problem_setup_finish_action)
1017  !Do nothing???
1018  CASE DEFAULT
1019  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1020  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1021  & " is invalid for a standard Stokes fluid."
1022  CALL flagerror(local_error,err,error,*999)
1023  END SELECT
1024  CASE(problem_setup_control_type)
1025  SELECT CASE(problem_setup%ACTION_TYPE)
1026  CASE(problem_setup_start_action)
1027  !Set up a simple control loop
1028  CALL control_loop_create_start(problem,control_loop,err,error,*999)
1029  CASE(problem_setup_finish_action)
1030  !Finish the control loops
1031  control_loop_root=>problem%CONTROL_LOOP
1032  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1033  CALL control_loop_create_finish(control_loop,err,error,*999)
1034  CASE DEFAULT
1035  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1036  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1037  & " is invalid for a standard Stokes fluid."
1038  CALL flagerror(local_error,err,error,*999)
1039  END SELECT
1040  CASE(problem_setup_solvers_type)
1041  !Get the control loop
1042  control_loop_root=>problem%CONTROL_LOOP
1043  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1044  SELECT CASE(problem_setup%ACTION_TYPE)
1045  CASE(problem_setup_start_action)
1046  !Start the solvers creation
1047  CALL solvers_create_start(control_loop,solvers,err,error,*999)
1048  CALL solvers_number_set(solvers,1,err,error,*999)
1049  !Set the solver to be a linear solver
1050  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1051  CALL solver_type_set(solver,solver_linear_type,err,error,*999)
1052  !Set solver defaults
1053  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
1054  CASE(problem_setup_finish_action)
1055  !Get the solvers
1056  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1057  !Finish the solvers creation
1058  CALL solvers_create_finish(solvers,err,error,*999)
1059  CASE DEFAULT
1060  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1061  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1062  & " is invalid for a standard Stokes fluid."
1063  CALL flagerror(local_error,err,error,*999)
1064  END SELECT
1065  CASE(problem_setup_solver_equations_type)
1066  SELECT CASE(problem_setup%ACTION_TYPE)
1067  CASE(problem_setup_start_action)
1068  !Get the control loop
1069  control_loop_root=>problem%CONTROL_LOOP
1070  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1071  !Get the solver
1072  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1073  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1074  !Create the solver equations
1075  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
1076  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
1077  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_static,err,error,*999)
1078  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
1079  CASE(problem_setup_finish_action)
1080  !Get the control loop
1081  control_loop_root=>problem%CONTROL_LOOP
1082  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1083  !Get the solver equations
1084  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1085  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1086  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1087  !Finish the solver equations creation
1088  CALL solver_equations_create_finish(solver_equations,err,error,*999)
1089  CASE DEFAULT
1090  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1091  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1092  & " is invalid for a standard Stokes fluid."
1093  CALL flagerror(local_error,err,error,*999)
1094  END SELECT
1095  CASE DEFAULT
1096  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1097  & " is invalid for a standard Stokes fluid."
1098  CALL flagerror(local_error,err,error,*999)
1099  END SELECT
1100  !Set problem subtype for transient Stokes problems
1101  CASE(problem_transient_stokes_subtype,problem_pgm_stokes_subtype)
1102  SELECT CASE(problem_setup%SETUP_TYPE)
1103  CASE(problem_setup_initial_type)
1104  SELECT CASE(problem_setup%ACTION_TYPE)
1105  CASE(problem_setup_start_action)
1106  !Do nothing????
1107  CASE(problem_setup_finish_action)
1108  !Do nothing????
1109  CASE DEFAULT
1110  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1111  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1112  & " is invalid for a transient Stokes fluid."
1113  CALL flagerror(local_error,err,error,*999)
1114  END SELECT
1115  CASE(problem_setup_control_type)
1116  SELECT CASE(problem_setup%ACTION_TYPE)
1117  CASE(problem_setup_start_action)
1118  !Set up a time control loop
1119  CALL control_loop_create_start(problem,control_loop,err,error,*999)
1120  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
1121  CASE(problem_setup_finish_action)
1122  !Finish the control loops
1123  control_loop_root=>problem%CONTROL_LOOP
1124  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1125  CALL control_loop_create_finish(control_loop,err,error,*999)
1126  CASE DEFAULT
1127  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1128  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1129  & " is invalid for a transient Stokes fluid."
1130  CALL flagerror(local_error,err,error,*999)
1131  END SELECT
1132  CASE(problem_setup_solvers_type)
1133  !Get the control loop
1134  control_loop_root=>problem%CONTROL_LOOP
1135  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1136  SELECT CASE(problem_setup%ACTION_TYPE)
1137  CASE(problem_setup_start_action)
1138  !Start the solvers creation
1139  CALL solvers_create_start(control_loop,solvers,err,error,*999)
1140  CALL solvers_number_set(solvers,1,err,error,*999)
1141  !Set the solver to be a first order dynamic solver
1142  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1143  CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
1144  CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
1145  !Set solver defaults
1146  CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
1147  CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
1148  CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
1149  CASE(problem_setup_finish_action)
1150  !Get the solvers
1151  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1152  !Finish the solvers creation
1153  CALL solvers_create_finish(solvers,err,error,*999)
1154  CASE DEFAULT
1155  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1156  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1157  & " is invalid for a transient Stokes fluid."
1158  CALL flagerror(local_error,err,error,*999)
1159  END SELECT
1160  CASE(problem_setup_solver_equations_type)
1161  SELECT CASE(problem_setup%ACTION_TYPE)
1162  CASE(problem_setup_start_action)
1163  !Get the control loop
1164  control_loop_root=>problem%CONTROL_LOOP
1165  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1166  !Get the solver
1167  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1168  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1169  !Create the solver equations
1170  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
1171  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
1172  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,&
1173  & err,error,*999)
1174  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
1175  CASE(problem_setup_finish_action)
1176  !Get the control loop
1177  control_loop_root=>problem%CONTROL_LOOP
1178  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1179  !Get the solver equations
1180  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1181  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1182  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1183  !Finish the solver equations creation
1184  CALL solver_equations_create_finish(solver_equations,err,error,*999)
1185  CASE DEFAULT
1186  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1187  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1188  & " is invalid for a transient Stokes fluid."
1189  CALL flagerror(local_error,err,error,*999)
1190  END SELECT
1191  CASE DEFAULT
1192  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1193  & " is invalid for a transient Stokes fluid."
1194  CALL flagerror(local_error,err,error,*999)
1195  END SELECT
1196  !Set problem subtype for ALE Stokes problems
1197  CASE(problem_ale_stokes_subtype)
1198  SELECT CASE(problem_setup%SETUP_TYPE)
1199  CASE(problem_setup_initial_type)
1200  SELECT CASE(problem_setup%ACTION_TYPE)
1201  CASE(problem_setup_start_action)
1202  !Do nothing????
1203  CASE(problem_setup_finish_action)
1204  !Do nothing????
1205  CASE DEFAULT
1206  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1207  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1208  & " is invalid for a ALE Stokes fluid."
1209  CALL flagerror(local_error,err,error,*999)
1210  END SELECT
1211  CASE(problem_setup_control_type)
1212  SELECT CASE(problem_setup%ACTION_TYPE)
1213  CASE(problem_setup_start_action)
1214  !Set up a time control loop
1215  CALL control_loop_create_start(problem,control_loop,err,error,*999)
1216  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
1217  CASE(problem_setup_finish_action)
1218  !Finish the control loops
1219  control_loop_root=>problem%CONTROL_LOOP
1220  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1221  CALL control_loop_create_finish(control_loop,err,error,*999)
1222  CASE DEFAULT
1223  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1224  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1225  & " is invalid for a ALE Stokes fluid."
1226  CALL flagerror(local_error,err,error,*999)
1227  END SELECT
1228  CASE(problem_setup_solvers_type)
1229  !Get the control loop
1230  control_loop_root=>problem%CONTROL_LOOP
1231  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1232  SELECT CASE(problem_setup%ACTION_TYPE)
1233  CASE(problem_setup_start_action)
1234  !Start the solvers creation
1235  CALL solvers_create_start(control_loop,solvers,err,error,*999)
1236  CALL solvers_number_set(solvers,2,err,error,*999)
1237  !Set the first solver to be a linear solver for the Laplace mesh movement problem
1238  CALL solvers_solver_get(solvers,1,mesh_solver,err,error,*999)
1239  CALL solver_type_set(mesh_solver,solver_linear_type,err,error,*999)
1240  !Set solver defaults
1241  CALL solver_library_type_set(mesh_solver,solver_petsc_library,err,error,*999)
1242  !Set the solver to be a first order dynamic solver
1243  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
1244  CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
1245  CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
1246  !Set solver defaults
1247  CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
1248  CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
1249  CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
1250  CASE(problem_setup_finish_action)
1251  !Get the solvers
1252  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1253  !Finish the solvers creation
1254  CALL solvers_create_finish(solvers,err,error,*999)
1255  CASE DEFAULT
1256  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1257  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1258  & " is invalid for a ALE Stokes fluid."
1259  CALL flagerror(local_error,err,error,*999)
1260  END SELECT
1261  CASE(problem_setup_solver_equations_type)
1262  SELECT CASE(problem_setup%ACTION_TYPE)
1263  CASE(problem_setup_start_action)
1264  !Get the control loop
1265  control_loop_root=>problem%CONTROL_LOOP
1266  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1267  !Get the solver
1268  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1269  CALL solvers_solver_get(solvers,1,mesh_solver,err,error,*999)
1270  !Create the solver equations
1271  CALL solver_equations_create_start(mesh_solver,mesh_solver_equations,err,error,*999)
1272  CALL solver_equations_linearity_type_set(mesh_solver_equations,solver_equations_linear,err,error,*999)
1273  CALL solver_equations_time_dependence_type_set(mesh_solver_equations,solver_equations_static,err,error,*999)
1274  CALL solver_equations_sparsity_type_set(mesh_solver_equations,solver_sparse_matrices,err,error,*999)
1275  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
1276  !Create the solver equations
1277  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
1278  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
1279  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,&
1280  & err,error,*999)
1281  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
1282  CASE(problem_setup_finish_action)
1283  !Get the control loop
1284  control_loop_root=>problem%CONTROL_LOOP
1285  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1286  !Get the solver equations
1287  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1288  CALL solvers_solver_get(solvers,1,mesh_solver,err,error,*999)
1289  CALL solver_solver_equations_get(mesh_solver,mesh_solver_equations,err,error,*999)
1290  !Finish the solver equations creation
1291  CALL solver_equations_create_finish(mesh_solver_equations,err,error,*999)
1292  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
1293  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1294  !Finish the solver equations creation
1295  CALL solver_equations_create_finish(solver_equations,err,error,*999)
1296  CASE DEFAULT
1297  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1298  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1299  & " is invalid for a ALE Stokes fluid."
1300  CALL flagerror(local_error,err,error,*999)
1301  END SELECT
1302  CASE DEFAULT
1303  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1304  & " is invalid for a ALE Stokes fluid."
1305  CALL flagerror(local_error,err,error,*999)
1306  END SELECT
1307  CASE DEFAULT
1308  local_error="Problem subtype "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
1309  & " is not valid for a Stokes fluid type of a fluid mechanics problem class."
1310  CALL flagerror(local_error,err,error,*999)
1311  END SELECT
1312  ELSE
1313  CALL flagerror("Problem is not associated.",err,error,*999)
1314  ENDIF
1315 
1316  exits("STOKES_PROBLEM_SETUP")
1317  RETURN
1318 999 errorsexits("STOKES_PROBLEM_SETUP",err,error)
1319  RETURN 1
1320  END SUBROUTINE stokes_problem_setup
1321 
1322 !
1323 !================================================================================================================================
1324 !
1325 
1327  SUBROUTINE stokes_finite_element_calculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
1328 
1329  !Argument variables
1330  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1331  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
1332  INTEGER(INTG), INTENT(OUT) :: ERR
1333  TYPE(varying_string), INTENT(OUT) :: ERROR
1334  !Local Variables
1335  INTEGER(INTG) FIELD_VAR_TYPE,ng,mh,mhs,mi,ms,nh,nhs,ni,ns,MESH_COMPONENT1,MESH_COMPONENT2, nhs_max, mhs_max, nhs_min, mhs_min
1336  REAL(DP) :: JGW,SUM,DXI_DX(3,3),PHIMS,PHINS,MU_PARAM,RHO_PARAM,DPHIMS_DXI(3),DPHINS_DXI(3)
1337  LOGICAL :: UPDATE_STIFFNESS_MATRIX, UPDATE_DAMPING_MATRIX,UPDATE_RHS_VECTOR
1338  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,DEPENDENT_BASIS1,DEPENDENT_BASIS2,GEOMETRIC_BASIS,INDEPENDENT_BASIS
1339  TYPE(equations_type), POINTER :: EQUATIONS
1340  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
1341  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
1342  TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
1343  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
1344  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
1345  TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
1346  TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
1347  TYPE(equations_matrix_type), POINTER :: STIFFNESS_MATRIX, DAMPING_MATRIX
1348  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD,INDEPENDENT_FIELD
1349  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
1350  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME,QUADRATURE_SCHEME1,QUADRATURE_SCHEME2
1351  TYPE(varying_string) :: LOCAL_ERROR
1352  INTEGER:: xv,out
1353  REAL(DP) :: AG_MATRIX(256,256) ! "A" Matrix ("G"radient part) - maximum size allocated
1354  REAL(DP) :: AL_MATRIX(256,256) ! "A" Matrix ("L"aplace part) - maximum size allocated
1355  REAL(DP) :: BT_MATRIX(256,256) ! "B" "T"ranspose Matrix - maximum size allocated
1356  REAL(DP) :: MT_MATRIX(256,256) ! "M"ass "T"ime Matrix - maximum size allocated
1357  REAL(DP) :: CT_MATRIX(256,256) ! "C"onvective "T"erm Matrix - maximum size allocated
1358  REAL(DP) :: ALE_MATRIX(256,256) ! "A"rbitrary "L"agrangian "E"ulerian Matrix - maximum size allocated
1359  REAL(DP) :: RH_VECTOR(256) ! "R"ight "H"and vector - maximum size allocated
1360  REAL(DP) :: W_VALUE(3)
1361  REAL(DP):: X(3)
1362 
1363 #ifdef TAUPROF
1364  CHARACTER(26) :: CVAR
1365  INTEGER :: GAUSS_POINT_LOOP_PHASE(2) = [ 0, 0 ]
1366  SAVE gauss_point_loop_phase
1367  INTEGER :: FIELD_COMPONENT_LOOP_PHASE(2) = [ 0, 0 ]
1368  SAVE field_component_loop_phase
1369 #endif
1370 
1371 !\todo: Reduce number of variables and parameters
1372 
1373  enters("STOKES_FINITE_ELEMENT_CALCULATE",err,error,*999)
1374 
1375  out=0
1376  ag_matrix=0.0_dp
1377  al_matrix=0.0_dp
1378  bt_matrix=0.0_dp
1379  mt_matrix=0.0_dp
1380  ct_matrix=0.0_dp
1381  ale_matrix=0.0_dp
1382  rh_vector=0.0_dp
1383  x=0.0_dp
1384 ! L=10.0_DP
1385 
1386  update_stiffness_matrix=.false.
1387  update_damping_matrix=.false.
1388  update_rhs_vector=.false.
1389 
1390  IF(ASSOCIATED(equations_set)) THEN
1391  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
1392  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
1393  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
1394  CALL flagerror("Equations set specification must have three entries for a Stokes type equations set.", &
1395  & err,error,*999)
1396  END IF
1397  equations=>equations_set%EQUATIONS
1398  IF(ASSOCIATED(equations)) THEN
1399  SELECT CASE(equations_set%SPECIFICATION(3))
1400  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype, &
1401  & equations_set_transient_stokes_subtype,problem_ale_stokes_subtype,problem_pgm_stokes_subtype)
1402  !Set Pointers
1403 #ifdef TAUPROF
1404  CALL tau_static_phase_start("SET POINTERS")
1405 #endif
1406  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
1407  independent_field=>equations%INTERPOLATION%INDEPENDENT_FIELD
1408  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
1409  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
1410  equations_matrices=>equations%EQUATIONS_MATRICES
1411  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
1412  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1413  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
1414  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1415  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
1416  rhs_vector=>equations_matrices%RHS_VECTOR
1417  equations_mapping=>equations%EQUATIONS_MAPPING
1418  SELECT CASE(equations_set%SPECIFICATION(3))
1419  CASE(equations_set_static_stokes_subtype,equations_set_laplace_stokes_subtype)
1420  linear_matrices=>equations_matrices%LINEAR_MATRICES
1421  stiffness_matrix=>linear_matrices%MATRICES(1)%PTR
1422  linear_mapping=>equations_mapping%LINEAR_MAPPING
1423  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1424  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1425  IF(ASSOCIATED(stiffness_matrix)) update_stiffness_matrix=stiffness_matrix%UPDATE_MATRIX
1426  IF(ASSOCIATED(rhs_vector)) update_rhs_vector=rhs_vector%UPDATE_VECTOR
1427  CASE(equations_set_transient_stokes_subtype)
1428  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
1429  stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
1430  damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
1431  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
1432  field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1433  field_var_type=field_variable%VARIABLE_TYPE
1434  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1435  damping_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1436  IF(ASSOCIATED(stiffness_matrix)) update_stiffness_matrix=stiffness_matrix%UPDATE_MATRIX
1437  IF(ASSOCIATED(damping_matrix)) update_damping_matrix=damping_matrix%UPDATE_MATRIX
1438  IF(ASSOCIATED(rhs_vector)) update_rhs_vector=rhs_vector%UPDATE_VECTOR
1439  CASE(equations_set_ale_stokes_subtype,equations_set_pgm_stokes_subtype)
1440  independent_field=>equations%INTERPOLATION%INDEPENDENT_FIELD
1441  independent_basis=>independent_field%DECOMPOSITION%DOMAIN(independent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)% &
1442  & ptr%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BASIS
1443  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
1444  stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
1445  damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
1446  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
1447  field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1448  field_var_type=field_variable%VARIABLE_TYPE
1449  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1450  damping_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1451  IF(ASSOCIATED(stiffness_matrix)) update_stiffness_matrix=stiffness_matrix%UPDATE_MATRIX
1452  IF(ASSOCIATED(damping_matrix)) update_damping_matrix=damping_matrix%UPDATE_MATRIX
1453  IF(ASSOCIATED(rhs_vector)) update_rhs_vector=rhs_vector%UPDATE_VECTOR
1454  CALL field_interpolation_parameters_element_get(field_mesh_velocity_set_type,element_number, &
1455  & equations%INTERPOLATION%INDEPENDENT_INTERP_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
1456  CASE DEFAULT
1457  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1458  & " is not valid for a Stokes fluid type of a fluid mechanics equations set class."
1459  CALL flagerror(local_error,err,error,*999)
1460  END SELECT
1461  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
1462  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
1463  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
1464  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
1465 #ifdef TAUPROF
1466  CALL tau_static_phase_stop("SET POINTERS")
1467 #endif
1468  !Start looping over Gauss points
1469  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
1470 #ifdef TAUPROF
1471  WRITE (cvar,'(a17,i2)') 'Gauss Point Loop ',ng
1472  CALL tau_phase_create_dynamic(gauss_point_loop_phase,cvar)
1473  CALL tau_phase_start(gauss_point_loop_phase)
1474  CALL tau_static_phase_start("INTERPOLATE")
1475 #endif
1476  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
1477  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
1478  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
1479  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
1480  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
1481  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
1482  IF(equations_set%SPECIFICATION(3)==equations_set_ale_stokes_subtype.OR. &
1483  & equations_set%SPECIFICATION(3)==equations_set_pgm_stokes_subtype) THEN
1484  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
1485  & independent_interp_point(field_u_variable_type)%PTR,err,error,*999)
1486  w_value(1)=equations%INTERPOLATION%INDEPENDENT_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
1487  w_value(2)=equations%INTERPOLATION%INDEPENDENT_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,no_part_deriv)
1488  IF(field_variable%NUMBER_OF_COMPONENTS==4) THEN
1489  w_value(3)=equations%INTERPOLATION%INDEPENDENT_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,no_part_deriv)
1490  END IF
1491  ELSE
1492  w_value=0.0_dp
1493  END IF
1494  !Define MU_PARAM, viscosity=1
1495  mu_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
1496  !Define RHO_PARAM, density=2
1497  rho_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,no_part_deriv)
1498 #ifdef TAUPROF
1499  CALL tau_static_phase_stop("INTERPOLATE")
1500 #endif
1501  !Calculate partial matrices
1502 !\todo: Check time spent here
1503  IF(equations_set%SPECIFICATION(3)==equations_set_static_stokes_subtype.OR. &
1504  & equations_set%SPECIFICATION(3)==equations_set_laplace_stokes_subtype.OR. &
1505  & equations_set%SPECIFICATION(3)==equations_set_ale_stokes_subtype.OR. &
1506  & equations_set%SPECIFICATION(3)==equations_set_pgm_stokes_subtype.OR. &
1507  & equations_set%SPECIFICATION(3)==equations_set_transient_stokes_subtype) THEN
1508  !Loop over field components
1509  mhs=0
1510  DO mh=1,(field_variable%NUMBER_OF_COMPONENTS-1)
1511 #ifdef TAUPROF
1512  WRITE (cvar,'(a22,i2)') 'Field Components Loop ',mh
1513  CALL tau_phase_create_dynamic(field_component_loop_phase,cvar)
1514  CALL tau_phase_start(field_component_loop_phase)
1515 
1516  CALL tau_static_phase_start("Field Set Variables")
1517 #endif
1518  mesh_component1=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1519  dependent_basis1=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component1)%PTR% &
1520  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1521  quadrature_scheme1=>dependent_basis1%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
1522  jgw=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
1523  & quadrature_scheme1%GAUSS_WEIGHTS(ng)
1524 #ifdef TAUPROF
1525  CALL tau_static_phase_stop("Field Set Variables")
1526  CALL tau_static_phase_start("Main Loop")
1527 #endif
1528  DO ms=1,dependent_basis1%NUMBER_OF_ELEMENT_PARAMETERS
1529  mhs=mhs+1
1530  nhs=0
1531  IF(update_stiffness_matrix.OR.update_damping_matrix) THEN
1532  !Loop over element columns
1533  DO nh=1,(field_variable%NUMBER_OF_COMPONENTS)
1534 !#ifdef TAUPROF
1535 ! CALL TAU_STATIC_PHASE_START("Element Set Variables")
1536 !#endif
1537  mesh_component2=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
1538  dependent_basis2=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component2)%PTR% &
1539  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1540  quadrature_scheme2=>dependent_basis2%QUADRATURE%QUADRATURE_SCHEME_MAP &
1541  & (basis_default_quadrature_scheme)%PTR
1542  ! JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME2%&
1543  ! &GAUSS_WEIGHTS(ng)
1544 !#ifdef TAUPROF
1545 ! CALL TAU_STATIC_PHASE_STOP("Element Set Variables")
1546 !#endif
1547  DO ns=1,dependent_basis2%NUMBER_OF_ELEMENT_PARAMETERS
1548  nhs=nhs+1
1549  !Calculate some variables used later on
1550 !\todo: Check time spent here
1551 !\todo: Re-introduce GU
1552 !#ifdef TAUPROF
1553 ! CALL TAU_STATIC_PHASE_START("PRECONDITIONS")
1554 !#endif
1555  DO ni=1,dependent_basis2%NUMBER_OF_XI
1556  DO mi=1,dependent_basis1%NUMBER_OF_XI
1557  dxi_dx(mi,ni)=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR% &
1558  & dxi_dx(mi,ni)
1559  END DO
1560  dphims_dxi(ni)=quadrature_scheme1%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
1561  dphins_dxi(ni)=quadrature_scheme2%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
1562  END DO !ni
1563  phims=quadrature_scheme1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
1564  phins=quadrature_scheme2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
1565  ! DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI
1566  ! DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI
1567  ! SUM=SUM-MU_PARAM*DPHIMSS_DXI(mi)*DPHINSS_DXI(ni)*EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%GU(mi,ni)
1568  ! ENDDO !ni
1569  ! ENDDO !mi
1570 !#ifdef TAUPROF
1571 ! CALL TAU_STATIC_PHASE_STOP("PRECONDITIONS")
1572 !#endif
1573  !Calculate Laplacian matrix
1574 !\todo: Check time spent here
1575 !\todo: Check IF statement
1576  IF(update_stiffness_matrix) THEN
1577 !#ifdef TAUPROF
1578 ! CALL TAU_STATIC_PHASE_START("A LAPLACE")
1579 !#endif
1580  !LAPLACE TYPE
1581  IF(nh==mh) THEN
1582  sum=0.0_dp
1583  !Calculate SUM
1584  DO xv=1,dependent_basis1%NUMBER_OF_XI
1585  DO mi=1,dependent_basis1%NUMBER_OF_XI
1586  DO ni=1,dependent_basis2%NUMBER_OF_XI
1587  sum=sum+mu_param*dphins_dxi(ni)*dxi_dx(ni,xv)*dphims_dxi(mi)*dxi_dx(mi,xv)
1588  ENDDO !ni
1589  ENDDO !mi
1590  ENDDO !x
1591  !Calculate MATRIX
1592  al_matrix(mhs,nhs)=al_matrix(mhs,nhs)+sum*jgw
1593  END IF
1594 !#ifdef TAUPROF
1595 ! CALL TAU_STATIC_PHASE_STOP("A LAPLACE")
1596 !#endif
1597  END IF
1598 !\todo: Check IF statement
1599  !Calculate standard matrix (gradient transpose type)
1600  IF(update_stiffness_matrix) THEN
1601 !#ifdef TAUPROF
1602 ! CALL TAU_STATIC_PHASE_START("A STANDARD")
1603 !#endif
1604  IF(equations_set%SPECIFICATION(3)/=equations_set_laplace_stokes_subtype) THEN
1605  IF(nh<field_variable%NUMBER_OF_COMPONENTS) THEN
1606  sum=0.0_dp
1607  !Calculate SUM
1608  DO mi=1,dependent_basis1%NUMBER_OF_XI
1609  DO ni=1,dependent_basis2%NUMBER_OF_XI
1610  !note mh/nh derivative in DXI_DX
1611  sum=sum+mu_param*dphins_dxi(mi)*dxi_dx(mi,mh)*dphims_dxi(ni)*dxi_dx(ni,nh)
1612  ENDDO !ni
1613  ENDDO !mi
1614  !Calculate MATRIX
1615  ag_matrix(mhs,nhs)=ag_matrix(mhs,nhs)+sum*jgw
1616  END IF
1617  END IF
1618 !#ifdef TAUPROF
1619 ! CALL TAU_STATIC_PHASE_STOP("A STANDARD")
1620 !#endif
1621  END IF
1622  !Calculate ALE matric contribution
1623 !\todo: Check time spent here
1624  IF(update_stiffness_matrix) THEN
1625 !#ifdef TAUPROF
1626 ! CALL TAU_STATIC_PHASE_START("ALE CONTRIBUTION")
1627 !#endif
1628  IF(equations_set%SPECIFICATION(3)==equations_set_ale_stokes_subtype.OR. &
1629  & equations_set%SPECIFICATION(3)==equations_set_pgm_stokes_subtype) THEN
1630  IF(nh==mh) THEN
1631  sum=0.0_dp
1632  !Calculate SUM
1633  DO mi=1,dependent_basis1%NUMBER_OF_XI
1634  DO ni=1,dependent_basis1%NUMBER_OF_XI
1635  sum=sum-rho_param*w_value(mi)*dphins_dxi(ni)*dxi_dx(ni,mi)*phims
1636  ENDDO !ni
1637  ENDDO !mi
1638  !Calculate MATRIX
1639  ale_matrix(mhs,nhs)=ale_matrix(mhs,nhs)+sum*jgw
1640  END IF
1641  END IF
1642 !#ifdef TAUPROF
1643 ! CALL TAU_STATIC_PHASE_STOP("ALE CONTRIBUTION")
1644 !#endif
1645  END IF
1646 !\todo: Check time spent here
1647  !Calculate pressure contribution (B transpose type)
1648  IF(update_stiffness_matrix) THEN
1649 !#ifdef TAUPROF
1650 ! CALL TAU_STATIC_PHASE_START("B TRANSPOSE")
1651 !#endif
1652  !LAPLACE TYPE
1653  IF(nh==field_variable%NUMBER_OF_COMPONENTS) THEN
1654  sum=0.0_dp
1655  !Calculate SUM
1656  DO ni=1,dependent_basis1%NUMBER_OF_XI
1657  sum=sum-phins*dphims_dxi(ni)*dxi_dx(ni,mh)
1658  ENDDO !ni
1659  !Calculate MATRIX
1660  bt_matrix(mhs,nhs)=bt_matrix(mhs,nhs)+sum*jgw
1661  END IF
1662 !#ifdef TAUPROF
1663 ! CALL TAU_STATIC_PHASE_STOP("B TRANSPOSE")
1664 !#endif
1665  END IF
1666  !Calculate mass matrix if needed
1667 !\todo: Check update flags
1668  IF(equations_set%SPECIFICATION(3)==equations_set_transient_stokes_subtype.OR. &
1669  & equations_set%SPECIFICATION(3)==equations_set_ale_stokes_subtype.OR. &
1670  & equations_set%SPECIFICATION(3)==equations_set_pgm_stokes_subtype) THEN
1671 !#ifdef TAUPROF
1672 ! CALL TAU_STATIC_PHASE_START("M")
1673 !#endif
1674  IF(update_damping_matrix) THEN
1675  IF(nh==mh) THEN
1676  sum=0.0_dp
1677  !Calculate SUM
1678  sum=phims*phins*rho_param
1679  !Calculate MATRIX
1680  mt_matrix(mhs,nhs)=mt_matrix(mhs,nhs)+sum*jgw
1681  END IF
1682  END IF
1683 !#ifdef TAUPROF
1684 ! CALL TAU_STATIC_PHASE_STOP("M")
1685 !#endif
1686  END IF
1687  ENDDO !ns
1688  ENDDO !nh
1689  ENDIF
1690  ENDDO !ms
1691 #ifdef TAUPROF
1692  CALL tau_static_phase_stop("Main Loop")
1693 
1694  CALL tau_phase_stop(field_component_loop_phase)
1695 #endif
1696  ENDDO !mh
1697 
1698 !#ifdef TAUPROF
1699 ! CALL TAU_STATIC_PHASE_STOP("FIELD COMPONENT LOOP")
1700 !#endif
1701 
1702  !Calculate analytic RHS
1703 #ifdef TAUPROF
1704  CALL tau_static_phase_start("RIGHT HAND SIDE FOR ANALYTIC SOLUTION")
1705 #endif
1706  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
1707  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_1.OR. &
1708  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_2.OR. &
1709  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_3.OR. &
1710  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_4.OR. &
1711  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_5.OR. &
1712  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_1.OR. &
1713  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_2.OR. &
1714  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_3.OR. &
1715  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_4.OR. &
1716  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_5) THEN
1717 
1718  mhs=0
1719  DO mh=1,(field_variable%NUMBER_OF_COMPONENTS-1)
1720  mesh_component1=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1721  dependent_basis1=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component1)%PTR% &
1722  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1723  quadrature_scheme1=>dependent_basis1%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
1724  jgw=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
1725  & quadrature_scheme1%GAUSS_WEIGHTS(ng)
1726  DO ms=1,dependent_basis1%NUMBER_OF_ELEMENT_PARAMETERS
1727  mhs=mhs+1
1728  phims=quadrature_scheme1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
1729  !note mh value derivative
1730  sum=0.0_dp
1731  x(1) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,1)
1732  x(2) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,1)
1733  IF(dependent_basis1%NUMBER_OF_XI==3) THEN
1734  x(3) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,1)
1735  END IF
1736  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_1) THEN
1737  IF(mh==1) THEN
1738  !Calculate SUM
1739  sum=0.0_dp
1740  ELSE IF(mh==2) THEN
1741  !Calculate SUM
1742  sum=phims*(-2.0_dp*mu_param/10.0_dp**2)
1743  ENDIF
1744  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_2) THEN
1745  IF(mh==1) THEN
1746  !Calculate SUM
1747  sum=0.0_dp
1748  ELSE IF(mh==2) THEN
1749  !Calculate SUM
1750  sum=phims*(-4.0_dp*mu_param/100.0_dp*exp((x(1)-x(2))/10.0_dp))
1751  ENDIF
1752  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_3) THEN
1753  IF(mh==1) THEN
1754  !Calculate SUM
1755  sum=0.0_dp
1756  ELSE IF(mh==2) THEN
1757  !Calculate SUM
1758  sum=phims*(16.0_dp*mu_param*pi*pi/100.0_dp*cos(2.0_dp*pi*x(2)/10.0_dp)*cos(2.0_dp*pi*x(1)/10.0_dp))
1759  ENDIF
1760  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_4) THEN
1761 ! do nothing!
1762  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_5) THEN
1763 ! do nothing!
1764  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_4) THEN
1765 ! do nothing!
1766  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_5) THEN
1767 ! do nothing!
1768  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_1) THEN
1769  IF(mh==1) THEN
1770  !Calculate SUM
1771  sum=0.0_dp
1772  ELSE IF(mh==2) THEN
1773  !Calculate SUM
1774  sum=phims*(-4.0_dp*mu_param/100.0_dp)
1775  ELSE IF(mh==3) THEN
1776  !Calculate SUM
1777  sum=phims*(-4.0_dp*mu_param/100.0_dp)
1778  ENDIF
1779  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_2) THEN
1780  IF(mh==1) THEN
1781  !Calculate SUM
1782  sum=0.0_dp
1783  ELSE IF(mh==2) THEN
1784  !Calculate SUM
1785  sum=phims*(-2.0_dp*mu_param/100.0_dp*(2.0_dp*exp((x(1)-x(2))/10.0_dp)+exp((x(2)-x(3))/10.0_dp)))
1786  ELSE IF(mh==3) THEN
1787  !Calculate SUM
1788  sum=phims*(-2.0_dp*mu_param/100.0_dp*(2.0_dp*exp((x(3)-x(1))/10.0_dp)+exp((x(2)-x(3))/10.0_dp)))
1789  ENDIF
1790  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_3) THEN
1791  IF(mh==1) THEN
1792  !Calculate SUM
1793  sum=0.0_dp
1794  ELSE IF(mh==2) THEN
1795  !Calculate SUM
1796  sum=phims*(36*mu_param*pi**2/100.0_dp*cos(2.0_dp*pi*x(2)/10.0_dp)*sin(2.0_dp*pi*x(3)/10.0_dp)* &
1797  & cos(2.0_dp*pi*x(1)/10.0_dp))
1798  ELSE IF(mh==3) THEN
1799  !Calculate SUM
1800  sum=0.0_dp
1801  ENDIF
1802  ENDIF
1803  !Calculate RH VECTOR
1804  rh_vector(mhs)=rh_vector(mhs)+sum*jgw
1805  ENDDO !ms
1806  ENDDO !mh
1807  ELSE
1808  rh_vector(mhs)=0.0_dp
1809  ENDIF
1810  ENDIF
1811 
1812 #ifdef TAUPROF
1813  CALL tau_static_phase_stop("RIGHT HAND SIDE FOR ANALYTIC SOLUTION")
1814 #endif
1815  END IF
1816 #ifdef TAUPROF
1817  CALL tau_phase_stop(gauss_point_loop_phase)
1818 #endif
1819  ENDDO !ng
1820  !Assemble matrices calculated above
1821 !\todo: Check time spent here
1822 #ifdef TAUPROF
1823  CALL tau_static_phase_start("ASSEMBLE STIFFNESS")
1824 #endif
1825  mhs_min=mhs
1826  mhs_max=nhs
1827  nhs_min=mhs
1828  nhs_max=nhs
1829  IF(equations_set%SPECIFICATION(3)==equations_set_static_stokes_subtype.OR. &
1830  & equations_set%SPECIFICATION(3)==equations_set_laplace_stokes_subtype.OR. &
1831  & equations_set%SPECIFICATION(3)==equations_set_ale_stokes_subtype.OR. &
1832  & equations_set%SPECIFICATION(3)==equations_set_pgm_stokes_subtype.OR. &
1833  & equations_set%SPECIFICATION(3)==equations_set_transient_stokes_subtype) THEN
1834  IF(update_stiffness_matrix) THEN
1835  stiffness_matrix%ELEMENT_MATRIX%MATRIX(1:mhs_min,1:nhs_min)=al_matrix(1:mhs_min,1:nhs_min)+ag_matrix(1:mhs_min, &
1836  & 1:nhs_min)+ale_matrix(1:mhs_min,1:nhs_min)
1837  stiffness_matrix%ELEMENT_MATRIX%MATRIX(1:mhs_min,nhs_min+1:nhs_max)=bt_matrix(1:mhs_min,nhs_min+1:nhs_max)
1838  DO mhs=mhs_min+1,mhs_max
1839  DO nhs=1,nhs_min
1840  !Transpose pressure type entries for mass equation
1841  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(nhs,mhs)
1842  END DO
1843  END DO
1844  END IF
1845  END IF
1846  IF(equations_set%SPECIFICATION(3)==equations_set_transient_stokes_subtype.OR. &
1847  & equations_set%SPECIFICATION(3)==equations_set_ale_stokes_subtype.OR. &
1848  & equations_set%SPECIFICATION(3)==equations_set_pgm_stokes_subtype) THEN
1849  IF(update_damping_matrix) THEN
1850  damping_matrix%ELEMENT_MATRIX%MATRIX(1:mhs_min,1:nhs_min)=mt_matrix(1:mhs_min,1:nhs_min)
1851  END IF
1852  END IF
1853 #ifdef TAUPROF
1854  CALL tau_static_phase_stop("ASSEMBLE STIFFNESS")
1855  CALL tau_static_phase_start("RIGHT HAND SIDE VECTOR")
1856 #endif
1857  !Assemble RHS vector
1858  IF(rhs_vector%FIRST_ASSEMBLY) THEN
1859  IF(update_rhs_vector) THEN
1860  rhs_vector%ELEMENT_VECTOR%VECTOR(1:mhs_max)=rh_vector(1:mhs_max)
1861  ENDIF
1862  ENDIF
1863 #ifdef TAUPROF
1864  CALL tau_static_phase_stop("RIGHT HAND SIDE VECTOR")
1865  CALL tau_static_phase_start("DEFINE SCALING")
1866 #endif
1867 !\todo: Check scaling factor for cubic hermite elements
1868  !Scale factor adjustment
1869  IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
1870  CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
1871  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
1872  mhs=0
1873  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
1874  !Loop over element rows
1875  mesh_component1=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1876  dependent_basis1=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component1)%PTR% &
1877  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1878  DO ms=1,dependent_basis1%NUMBER_OF_ELEMENT_PARAMETERS
1879  mhs=mhs+1
1880  nhs=0
1881  IF(update_stiffness_matrix.OR.update_damping_matrix) THEN
1882  !Loop over element columns
1883  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
1884  mesh_component2=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
1885  dependent_basis2=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component2)%PTR% &
1886  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1887  DO ns=1,dependent_basis2%NUMBER_OF_ELEMENT_PARAMETERS
1888  nhs=nhs+1
1889  IF(update_stiffness_matrix)THEN
1890  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
1891  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
1892  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
1893  END IF
1894  IF(update_damping_matrix)THEN
1895  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
1896  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
1897  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
1898  END IF
1899  ENDDO !ns
1900  ENDDO !nh
1901  ENDIF
1902  IF(update_rhs_vector) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
1903  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
1904  ENDDO !ms
1905  ENDDO !mh
1906  ENDIF
1907 #ifdef TAUPROF
1908  CALL tau_static_phase_stop("DEFINE SCALING")
1909 #endif
1910  CASE DEFAULT
1911  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1912  & " is not valid for a Stokes fluid type of a fluid mechanics equations set class."
1913  CALL flagerror(local_error,err,error,*999)
1914  END SELECT
1915  ELSE
1916  CALL flagerror("Equations set equations is not associated.",err,error,*999)
1917  ENDIF
1918  ELSE
1919  CALL flagerror("Equations set is not associated.",err,error,*999)
1920  ENDIF
1921 
1922  exits("STOKES_FINITE_ELEMENT_CALCULATE")
1923  RETURN
1924 999 errorsexits("STOKES_FINITE_ELEMENT_CALCULATE",err,error)
1925  RETURN 1
1926  END SUBROUTINE stokes_finite_element_calculate
1927 
1928  !
1929  !================================================================================================================================
1930  !
1931 
1933  SUBROUTINE stokes_post_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
1934 
1935  !Argument variables
1936  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
1937  TYPE(solver_type), POINTER :: SOLVER
1938  INTEGER(INTG), INTENT(OUT) :: ERR
1939  TYPE(varying_string), INTENT(OUT) :: ERROR
1940  !Local Variables
1941  TYPE(solver_type), POINTER :: SOLVER2
1942  TYPE(varying_string) :: LOCAL_ERROR
1943 
1944  enters("STOKES_POST_SOLVE",err,error,*999)
1945  NULLIFY(solver2)
1946 
1947  IF(ASSOCIATED(control_loop)) THEN
1948  IF(ASSOCIATED(solver)) THEN
1949  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
1950  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
1951  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1952  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
1953  CALL flagerror("Problem specification must have three entries for a Stokes problem.",err,error,*999)
1954  END IF
1955  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
1956  CASE(problem_static_stokes_subtype,problem_laplace_stokes_subtype)
1957  CALL stokes_post_solve_output_data(control_loop,solver,err,error,*999)
1958  CASE(problem_pgm_stokes_subtype)
1959  CALL stokes_post_solve_output_data(control_loop,solver,err,error,*999)
1960  CASE(problem_transient_stokes_subtype)
1961  CALL stokes_post_solve_output_data(control_loop,solver,err,error,*999)
1962  CASE(problem_ale_stokes_subtype)
1963  !Post solve for the linear solver
1964  IF(solver%SOLVE_TYPE==solver_linear_type) THEN
1965  CALL write_string(general_output_type,"Mesh movement post solve... ",err,error,*999)
1966  CALL solvers_solver_get(solver%SOLVERS,2,solver2,err,error,*999)
1967  IF(ASSOCIATED(solver2%DYNAMIC_SOLVER)) THEN
1968  solver2%DYNAMIC_SOLVER%ALE=.true.
1969  ELSE
1970  CALL flagerror("Dynamic solver is not associated for ALE problem.",err,error,*999)
1971  END IF
1972  !Post solve for the linear solver
1973  ELSE IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
1974  CALL write_string(general_output_type,"ALE Stokes post solve... ",err,error,*999)
1975  CALL stokes_post_solve_output_data(control_loop,solver,err,error,*999)
1976  END IF
1977  CASE DEFAULT
1978  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
1979  & " is not valid for a Stokes fluid type of a fluid mechanics problem class."
1980  CALL flagerror(local_error,err,error,*999)
1981  END SELECT
1982  ELSE
1983  CALL flagerror("Problem is not associated.",err,error,*999)
1984  ENDIF
1985  ELSE
1986  CALL flagerror("Solver is not associated.",err,error,*999)
1987  ENDIF
1988  ELSE
1989  CALL flagerror("Control loop is not associated.",err,error,*999)
1990  ENDIF
1991 
1992  exits("STOKES_POST_SOLVE")
1993  RETURN
1994 999 errorsexits("STOKES_POST_SOLVE",err,error)
1995  RETURN 1
1996  END SUBROUTINE stokes_post_solve
1997 
1998  !
1999  !================================================================================================================================
2000  !
2001 
2003  SUBROUTINE stokes_pre_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
2004 
2005  !Argument variables
2006  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
2007  TYPE(solver_type), POINTER :: SOLVER
2008  INTEGER(INTG), INTENT(OUT) :: ERR
2009  TYPE(varying_string), INTENT(OUT) :: ERROR
2010  !Local Variables
2011  TYPE(solver_type), POINTER :: SOLVER2
2012  TYPE(varying_string) :: LOCAL_ERROR
2013 
2014  enters("STOKES_PRE_SOLVE",err,error,*999)
2015  NULLIFY(solver2)
2016 
2017  IF(ASSOCIATED(control_loop)) THEN
2018  IF(ASSOCIATED(solver)) THEN
2019  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
2020  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
2021  CALL flagerror("Problem specification is not allocated.",err,error,*999)
2022  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
2023  CALL flagerror("Problem specification must have three entries for a Stokes problem.",err,error,*999)
2024  END IF
2025  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
2026  CASE(problem_static_stokes_subtype,problem_laplace_stokes_subtype)
2027  ! do nothing ???
2028  CASE(problem_transient_stokes_subtype)
2029  ! do nothing ???
2030  CALL stokes_pre_solve_update_boundary_conditions(control_loop,solver,err,error,*999)
2031  CASE(problem_pgm_stokes_subtype)
2032  ! do nothing ???
2033  !First update mesh and calculates boundary velocity values
2034  CALL stokes_pre_solve_ale_update_mesh(control_loop,solver,err,error,*999)
2035  !Then apply both normal and moving mesh boundary conditions
2036  CALL stokes_pre_solve_update_boundary_conditions(control_loop,solver,err,error,*999)
2037  CASE(problem_ale_stokes_subtype)
2038  !Pre solve for the linear solver
2039  IF(solver%SOLVE_TYPE==solver_linear_type) THEN
2040  CALL write_string(general_output_type,"Mesh movement pre solve... ",err,error,*999)
2041  !Update boundary conditions for mesh-movement
2042  CALL stokes_pre_solve_update_boundary_conditions(control_loop,solver,err,error,*999)
2043  CALL solvers_solver_get(solver%SOLVERS,2,solver2,err,error,*999)
2044  IF(ASSOCIATED(solver2%DYNAMIC_SOLVER)) THEN
2045 !\todo: Avoid ALE flag in future
2046  solver2%DYNAMIC_SOLVER%ALE=.false.
2047  ELSE
2048  CALL flagerror("Dynamic solver is not associated for ALE problem.",err,error,*999)
2049  END IF
2050  !Update material properties for Laplace mesh movement
2051  CALL stokes_pre_solve_ale_update_parameters(control_loop,solver,err,error,*999)
2052  !Pre solve for the linear solver
2053  ELSE IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
2054  CALL write_string(general_output_type,"ALE Stokes pre solve... ",err,error,*999)
2055  IF(solver%DYNAMIC_SOLVER%ALE) THEN
2056  !First update mesh and calculates boundary velocity values
2057  CALL stokes_pre_solve_ale_update_mesh(control_loop,solver,err,error,*999)
2058  !Then apply both normal and moving mesh boundary conditions
2059  CALL stokes_pre_solve_update_boundary_conditions(control_loop,solver,err,error,*999)
2060  ELSE
2061  CALL flagerror("Mesh motion calculation not successful for ALE problem.",err,error,*999)
2062  END IF
2063  ELSE
2064  CALL flagerror("Solver type is not associated for ALE problem.",err,error,*999)
2065  END IF
2066  CASE DEFAULT
2067  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
2068  & " is not valid for a Stokes fluid type of a fluid mechanics problem class."
2069  CALL flagerror(local_error,err,error,*999)
2070  END SELECT
2071  ELSE
2072  CALL flagerror("Problem is not associated.",err,error,*999)
2073  ENDIF
2074  ELSE
2075  CALL flagerror("Solver is not associated.",err,error,*999)
2076  ENDIF
2077  ELSE
2078  CALL flagerror("Control loop is not associated.",err,error,*999)
2079  ENDIF
2080 
2081  exits("STOKES_PRE_SOLVE")
2082  RETURN
2083 999 errorsexits("STOKES_PRE_SOLVE",err,error)
2084  RETURN 1
2085  END SUBROUTINE stokes_pre_solve
2086  !
2087  !================================================================================================================================
2088  !
2089 
2091  SUBROUTINE stokes_pre_solve_update_boundary_conditions(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
2092 
2093  !Argument variables
2094  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
2095  TYPE(solver_type), POINTER :: SOLVER
2096  INTEGER(INTG), INTENT(OUT) :: ERR
2097  TYPE(varying_string), INTENT(OUT) :: ERROR
2098  !Local Variables
2099  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
2100  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
2101  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2102  TYPE(equations_type), POINTER :: EQUATIONS
2103  TYPE(varying_string) :: LOCAL_ERROR
2104  TYPE(boundary_conditions_variable_type), POINTER :: BOUNDARY_CONDITIONS_VARIABLE
2105  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
2106  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
2107  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
2108  TYPE(domain_type), POINTER :: DOMAIN
2109  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
2110  TYPE(field_interpolated_point_ptr_type), POINTER :: INTERPOLATED_POINT(:)
2111  TYPE(field_interpolation_parameters_ptr_type), POINTER :: INTERPOLATION_PARAMETERS(:)
2112  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT,DISPLACEMENT_VALUE,VALUE,XI_COORDINATES(3)
2113  REAL(DP) :: T_COORDINATES(20,3)
2114  INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,BOUNDARY_CONDITION_CHECK_VARIABLE,GLOBAL_DERIV_INDEX,node_idx,variable_type
2115  INTEGER(INTG) :: variable_idx,local_ny,ANALYTIC_FUNCTION_TYPE,component_idx,deriv_idx,dim_idx
2116  INTEGER(INTG) :: element_idx,en_idx,I,J,K,number_of_nodes_xic(3)
2117  REAL(DP) :: X(3),MU_PARAM,RHO_PARAM
2118  REAL(DP), POINTER :: MESH_VELOCITY_VALUES(:), GEOMETRIC_PARAMETERS(:)
2119  REAL(DP), POINTER :: BOUNDARY_VALUES(:)
2120 
2121 
2122  enters("STOKES_PRE_SOLVE_UPDATE_BOUNDARY_CONDITIONS",err,error,*999)
2123 
2124  IF(ASSOCIATED(control_loop)) THEN
2125  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
2126 ! write(*,*)'CURRENT_TIME = ',CURRENT_TIME
2127 ! write(*,*)'TIME_INCREMENT = ',TIME_INCREMENT
2128  IF(ASSOCIATED(solver)) THEN
2129  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
2130  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
2131  CALL flagerror("Problem specification is not allocated.",err,error,*999)
2132  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
2133  CALL flagerror("Problem specification must have three entries for a Stokes problem.",err,error,*999)
2134  END IF
2135  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
2136  CASE(problem_static_stokes_subtype,problem_laplace_stokes_subtype)
2137  ! do nothing ???
2138  CASE(problem_transient_stokes_subtype)
2139  solver_equations=>solver%SOLVER_EQUATIONS
2140  IF(ASSOCIATED(solver_equations)) THEN
2141  solver_mapping=>solver_equations%SOLVER_MAPPING
2142  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
2143  IF(ASSOCIATED(equations)) THEN
2144  equations_set=>equations%EQUATIONS_SET
2145  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
2146  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_4 .OR. &
2147  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_1 .OR. &
2148  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_4 .OR. &
2149  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_three_dim_5 .OR. &
2150  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_stokes_equation_two_dim_5) THEN
2151  IF(ASSOCIATED(equations_set)) THEN
2152  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
2153  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
2154  IF(ASSOCIATED(dependent_field)) THEN
2155  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
2156  IF(ASSOCIATED(geometric_field)) THEN
2157  CALL field_number_of_components_get(geometric_field,field_u_variable_type,&
2158  & number_of_dimensions,err,error,*999)
2159  NULLIFY(interpolation_parameters)
2160  NULLIFY(interpolated_point)
2161  CALL field_interpolation_parameters_initialise(geometric_field,interpolation_parameters,err,error, &
2162  & *999)
2163  CALL field_interpolated_points_initialise(interpolation_parameters,interpolated_point,err,error,*999)
2164  NULLIFY(geometric_variable)
2165  CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
2166  NULLIFY(geometric_parameters)
2167  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,&
2168  & geometric_parameters,err,error,*999)
2169  DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES
2170  variable_type=dependent_field%VARIABLES(variable_idx)%VARIABLE_TYPE
2171  field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
2172  IF(ASSOCIATED(field_variable)) THEN
2173  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
2174  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE== &
2175  & field_node_based_interpolation) THEN
2176  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
2177  IF(ASSOCIATED(domain)) THEN
2178  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
2179  domain_nodes=>domain%TOPOLOGY%NODES
2180  IF(ASSOCIATED(domain_nodes)) THEN
2181  !Should be replaced by boundary node flag
2182  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
2183  element_idx=domain%topology%nodes%nodes(node_idx)%surrounding_elements(1)
2184  CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
2185  & interpolation_parameters(field_u_variable_type)%PTR,err,error,*999)
2186  en_idx=0
2187  xi_coordinates=0.0_dp
2188  number_of_nodes_xic(1)=domain%topology%elements%elements(element_idx)% &
2189  & basis%number_of_nodes_xic(1)
2190  number_of_nodes_xic(2)=domain%topology%elements%elements(element_idx)% &
2191  & basis%number_of_nodes_xic(2)
2192  IF(number_of_dimensions==3) THEN
2193  number_of_nodes_xic(3)=domain%topology%elements%elements(element_idx)%basis% &
2194  & number_of_nodes_xic(3)
2195  ELSE
2196  number_of_nodes_xic(3)=1
2197  ENDIF
2198 !\todo: change definitions as soon as adjacent elements / boundary elements calculation works for simplex
2199  IF(domain%topology%elements%maximum_number_of_element_parameters==4.OR. &
2200  & domain%topology%elements%maximum_number_of_element_parameters==9.OR. &
2201  & domain%topology%elements%maximum_number_of_element_parameters==16.OR. &
2202  & domain%topology%elements%maximum_number_of_element_parameters==8.OR. &
2203  & domain%topology%elements%maximum_number_of_element_parameters==27.OR. &
2204  & domain%topology%elements%maximum_number_of_element_parameters==64) THEN
2205  DO k=1,number_of_nodes_xic(3)
2206  DO j=1,number_of_nodes_xic(2)
2207  DO i=1,number_of_nodes_xic(1)
2208  en_idx=en_idx+1
2209  IF(domain%topology%elements%elements(element_idx)% &
2210  & element_nodes(en_idx)==node_idx) EXIT
2211  xi_coordinates(1)=xi_coordinates(1)+(1.0_dp/(number_of_nodes_xic(1)-1))
2212  ENDDO
2213  IF(domain%topology%elements%elements(element_idx)% &
2214  & element_nodes(en_idx)==node_idx) EXIT
2215  xi_coordinates(1)=0.0_dp
2216  xi_coordinates(2)=xi_coordinates(2)+(1.0_dp/(number_of_nodes_xic(2)-1))
2217  ENDDO
2218  IF(domain%topology%elements%elements(element_idx)% &
2219  & element_nodes(en_idx)==node_idx) EXIT
2220  xi_coordinates(1)=0.0_dp
2221  xi_coordinates(2)=0.0_dp
2222  IF(number_of_nodes_xic(3)/=1) THEN
2223  xi_coordinates(3)=xi_coordinates(3)+(1.0_dp/(number_of_nodes_xic(3)-1))
2224  ENDIF
2225  ENDDO
2226  CALL field_interpolate_xi(no_part_deriv,xi_coordinates, &
2227  & interpolated_point(field_u_variable_type)%PTR,err,error,*999)
2228  ELSE
2229 !\todo: Use boundary flag
2230  IF(domain%topology%elements%maximum_number_of_element_parameters==3) THEN
2231  t_coordinates(1,1:2)=[0.0_dp,1.0_dp]
2232  t_coordinates(2,1:2)=[1.0_dp,0.0_dp]
2233  t_coordinates(3,1:2)=[1.0_dp,1.0_dp]
2234  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==6) THEN
2235  t_coordinates(1,1:2)=[0.0_dp,1.0_dp]
2236  t_coordinates(2,1:2)=[1.0_dp,0.0_dp]
2237  t_coordinates(3,1:2)=[1.0_dp,1.0_dp]
2238  t_coordinates(4,1:2)=[0.5_dp,0.5_dp]
2239  t_coordinates(5,1:2)=[1.0_dp,0.5_dp]
2240  t_coordinates(6,1:2)=[0.5_dp,1.0_dp]
2241  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==10.AND. &
2242  & number_of_dimensions==2) THEN
2243  t_coordinates(1,1:2)=[0.0_dp,1.0_dp]
2244  t_coordinates(2,1:2)=[1.0_dp,0.0_dp]
2245  t_coordinates(3,1:2)=[1.0_dp,1.0_dp]
2246  t_coordinates(4,1:2)=[1.0_dp/3.0_dp,2.0_dp/3.0_dp]
2247  t_coordinates(5,1:2)=[2.0_dp/3.0_dp,1.0_dp/3.0_dp]
2248  t_coordinates(6,1:2)=[1.0_dp,1.0_dp/3.0_dp]
2249  t_coordinates(7,1:2)=[1.0_dp,2.0_dp/3.0_dp]
2250  t_coordinates(8,1:2)=[2.0_dp/3.0_dp,1.0_dp]
2251  t_coordinates(9,1:2)=[1.0_dp/3.0_dp,1.0_dp]
2252  t_coordinates(10,1:2)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp]
2253  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==4) THEN
2254  t_coordinates(1,1:3)=[0.0_dp,1.0_dp,1.0_dp]
2255  t_coordinates(2,1:3)=[1.0_dp,0.0_dp,1.0_dp]
2256  t_coordinates(3,1:3)=[1.0_dp,1.0_dp,0.0_dp]
2257  t_coordinates(4,1:3)=[1.0_dp,1.0_dp,1.0_dp]
2258  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==10.AND. &
2259  & number_of_dimensions==3) THEN
2260  t_coordinates(1,1:3)=[0.0_dp,1.0_dp,1.0_dp]
2261  t_coordinates(2,1:3)=[1.0_dp,0.0_dp,1.0_dp]
2262  t_coordinates(3,1:3)=[1.0_dp,1.0_dp,0.0_dp]
2263  t_coordinates(4,1:3)=[1.0_dp,1.0_dp,1.0_dp]
2264  t_coordinates(5,1:3)=[0.5_dp,0.5_dp,1.0_dp]
2265  t_coordinates(6,1:3)=[0.5_dp,1.0_dp,0.5_dp]
2266  t_coordinates(7,1:3)=[0.5_dp,1.0_dp,1.0_dp]
2267  t_coordinates(8,1:3)=[1.0_dp,0.5_dp,0.5_dp]
2268  t_coordinates(9,1:3)=[1.0_dp,1.0_dp,0.5_dp]
2269  t_coordinates(10,1:3)=[1.0_dp,0.5_dp,1.0_dp]
2270  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==20) THEN
2271  t_coordinates(1,1:3)=[0.0_dp,1.0_dp,1.0_dp]
2272  t_coordinates(2,1:3)=[1.0_dp,0.0_dp,1.0_dp]
2273  t_coordinates(3,1:3)=[1.0_dp,1.0_dp,0.0_dp]
2274  t_coordinates(4,1:3)=[1.0_dp,1.0_dp,1.0_dp]
2275  t_coordinates(5,1:3)=[1.0_dp/3.0_dp,2.0_dp/3.0_dp,1.0_dp]
2276  t_coordinates(6,1:3)=[2.0_dp/3.0_dp,1.0_dp/3.0_dp,1.0_dp]
2277  t_coordinates(7,1:3)=[1.0_dp/3.0_dp,1.0_dp,2.0_dp/3.0_dp]
2278  t_coordinates(8,1:3)=[2.0_dp/3.0_dp,1.0_dp,1.0_dp/3.0_dp]
2279  t_coordinates(9,1:3)=[1.0_dp/3.0_dp,1.0_dp,1.0_dp]
2280  t_coordinates(10,1:3)=[2.0_dp/3.0_dp,1.0_dp,1.0_dp]
2281  t_coordinates(11,1:3)=[1.0_dp,1.0_dp/3.0_dp,2.0_dp/3.0_dp]
2282  t_coordinates(12,1:3)=[1.0_dp,2.0_dp/3.0_dp,1.0_dp/3.0_dp]
2283  t_coordinates(13,1:3)=[1.0_dp,1.0_dp,1.0_dp/3.0_dp]
2284  t_coordinates(14,1:3)=[1.0_dp,1.0_dp,2.0_dp/3.0_dp]
2285  t_coordinates(15,1:3)=[1.0_dp,1.0_dp/3.0_dp,1.0_dp]
2286  t_coordinates(16,1:3)=[1.0_dp,2.0_dp/3.0_dp,1.0_dp]
2287  t_coordinates(17,1:3)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp,2.0_dp/3.0_dp]
2288  t_coordinates(18,1:3)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp,1.0_dp]
2289  t_coordinates(19,1:3)=[2.0_dp/3.0_dp,1.0_dp,2.0_dp/3.0_dp]
2290  t_coordinates(20,1:3)=[1.0_dp,2.0_dp/3.0_dp,2.0_dp/3.0_dp]
2291  ENDIF
2292  DO k=1,domain%topology%elements%maximum_number_of_element_parameters
2293  IF(domain%topology%elements%elements(element_idx)%element_nodes(k)==node_idx) EXIT
2294  ENDDO
2295  IF(number_of_dimensions==2) THEN
2296  CALL field_interpolate_xi(no_part_deriv,t_coordinates(k,1:2), &
2297  & interpolated_point(field_u_variable_type)%PTR,err,error,*999)
2298  ELSE IF(number_of_dimensions==3) THEN
2299  CALL field_interpolate_xi(no_part_deriv,t_coordinates(k,1:3), &
2300  & interpolated_point(field_u_variable_type)%PTR,err,error,*999)
2301  ENDIF
2302  ENDIF
2303  x=0.0_dp
2304  DO dim_idx=1,number_of_dimensions
2305  x(dim_idx)=interpolated_point(field_u_variable_type)%PTR%VALUES(dim_idx,1)
2306  ENDDO !dim_idx
2307  !Loop over the derivatives
2308  CALL boundary_conditions_variable_get(solver_equations%BOUNDARY_CONDITIONS, &
2309  & dependent_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR, &
2310  & boundary_conditions_variable,err,error,*999)
2311  IF(ASSOCIATED(boundary_conditions_variable)) THEN
2312  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
2313  analytic_function_type=equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE
2314  global_deriv_index=domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)% &
2315  & global_derivative_index
2316  materials_field=>equations_set%MATERIALS%MATERIALS_FIELD
2317  !Define MU_PARAM, density=1
2318  mu_param=materials_field%variables(1)%parameter_sets%parameter_sets(1)%ptr% &
2319  & parameters%cmiss%data_dp(1)
2320  !Define RHO_PARAM, density=2
2321  IF(analytic_function_type==equations_set_stokes_equation_two_dim_4.OR. &
2322  & analytic_function_type==equations_set_stokes_equation_two_dim_5.OR. &
2323  & analytic_function_type==equations_set_stokes_equation_three_dim_4.OR. &
2324  & analytic_function_type==equations_set_stokes_equation_three_dim_5) THEN
2325  rho_param=materials_field%variables(1)%parameter_sets%parameter_sets(1)%ptr% &
2326  & parameters%cmiss%data_dp(2)
2327  ELSE
2328  rho_param=0.0_dp
2329  ENDIF
2330  CALL stokes_equation_analytic_functions(VALUE,x,mu_param,rho_param,current_time, &
2331  & variable_type, &
2332  & global_deriv_index,analytic_function_type,number_of_dimensions, &
2333  & field_variable%NUMBER_OF_COMPONENTS,component_idx,err,error,*999)
2334  !Default to version 1 of each node derivative
2335  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
2336  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
2337  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
2338  & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
2339  boundary_condition_check_variable=boundary_conditions_variable% &
2340  & condition_types(local_ny)
2341  IF(boundary_condition_check_variable==boundary_condition_fixed) THEN
2342  CALL field_parameter_set_update_local_dof(dependent_field, &
2343  & variable_type,field_values_set_type,local_ny, &
2344  & VALUE,err,error,*999)
2345  ENDIF
2346  ENDDO !deriv_idx
2347  ELSE
2348  CALL flagerror("Boundary conditions U variable is not associated",err,error,*999)
2349  ENDIF
2350  ENDDO !node_idx
2351  ELSE
2352  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
2353  ENDIF
2354  ELSE
2355  CALL flagerror("Domain topology is not associated.",err,error,*999)
2356  ENDIF
2357  ELSE
2358  CALL flagerror("Domain is not associated.",err,error,*999)
2359  ENDIF
2360  ELSE
2361  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
2362  ENDIF
2363  ENDDO !component_idx
2364  CALL field_parameter_set_update_start(dependent_field,variable_type, &
2365  & field_analytic_values_set_type,err,error,*999)
2366  CALL field_parameter_set_update_finish(dependent_field,variable_type, &
2367  & field_analytic_values_set_type,err,error,*999)
2368  CALL field_parameter_set_update_start(dependent_field,variable_type, &
2369  & field_values_set_type,err,error,*999)
2370  CALL field_parameter_set_update_finish(dependent_field,variable_type, &
2371  & field_values_set_type,err,error,*999)
2372  ELSE
2373  CALL flagerror("Field variable is not associated.",err,error,*999)
2374  ENDIF
2375  ENDDO !variable_idx
2376  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,&
2377  & field_values_set_type,geometric_parameters,err,error,*999)
2378  ELSE
2379  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
2380  ENDIF
2381  ELSE
2382  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
2383  ENDIF
2384  ELSE
2385  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
2386  ENDIF
2387  ELSE
2388  CALL flagerror("Equations set is not associated.",err,error,*999)
2389  ENDIF
2390  ENDIF
2391  ENDIF
2392  ELSE
2393  CALL flagerror("Equations are not associated.",err,error,*999)
2394  END IF
2395  ELSE
2396  CALL flagerror("Solver equations are not associated.",err,error,*999)
2397  END IF
2398  CALL field_parameter_set_update_start(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2399  & field_values_set_type,err,error,*999)
2400  CALL field_parameter_set_update_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2401  & field_values_set_type,err,error,*999)
2402  CASE(problem_pgm_stokes_subtype)
2403  !Pre solve for the dynamic solver
2404  IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
2405  CALL write_string(general_output_type,"Mesh movement change boundary conditions... ",err,error,*999)
2406  solver_equations=>solver%SOLVER_EQUATIONS
2407  IF(ASSOCIATED(solver_equations)) THEN
2408  solver_mapping=>solver_equations%SOLVER_MAPPING
2409  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
2410  IF(ASSOCIATED(equations)) THEN
2411  equations_set=>equations%EQUATIONS_SET
2412  IF(ASSOCIATED(equations_set)) THEN
2413  boundary_conditions=>solver_equations%BOUNDARY_CONDITIONS
2414  IF(ASSOCIATED(boundary_conditions)) THEN
2415  CALL boundary_conditions_variable_get(boundary_conditions,equations_set%DEPENDENT%DEPENDENT_FIELD% &
2416  & variable_type_map(field_u_variable_type)%PTR,boundary_conditions_variable,err,error,*999)
2417  IF(ASSOCIATED(boundary_conditions_variable)) THEN
2418  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2419  & number_of_dimensions,err,error,*999)
2420  NULLIFY(mesh_velocity_values)
2421  CALL field_parameter_set_data_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
2422  & field_mesh_velocity_set_type,mesh_velocity_values,err,error,*999)
2423  NULLIFY(boundary_values)
2424  CALL field_parameter_set_data_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
2425  & field_boundary_set_type,boundary_values,err,error,*999)
2426  CALL fluid_mechanics_io_read_boundary_conditions(solver_linear_type,boundary_values, &
2427  & number_of_dimensions,boundary_condition_fixed_inlet,control_loop%TIME_LOOP%INPUT_NUMBER, &
2428  & control_loop%TIME_LOOP%ITERATION_NUMBER,current_time,1.0_dp)
2429 ! DO equations_row_number=1,EQUATIONS%EQUATIONS_MAPPING%TOTAL_NUMBER_OF_ROWS
2430 ! xxxxxxxxxxxxxxxxxxxxxx
2431  DO variable_idx=1,equations_set%DEPENDENT%DEPENDENT_FIELD%NUMBER_OF_VARIABLES
2432  variable_type=equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
2433  field_variable=>equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
2434  IF(ASSOCIATED(field_variable)) THEN
2435  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
2436  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
2437  IF(ASSOCIATED(domain)) THEN
2438  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
2439  domain_nodes=>domain%TOPOLOGY%NODES
2440  IF(ASSOCIATED(domain_nodes)) THEN
2441  !Loop over the local nodes excluding the ghosts.
2442  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
2443  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
2444  !Default to version 1 of each node derivative
2445  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
2446  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
2447 ! xxxxxxxxxxxxxxxxxxxxxxxxx
2448  displacement_value=0.0_dp
2449  boundary_condition_check_variable=boundary_conditions_variable% &
2450  & condition_types(local_ny)
2451  IF(boundary_condition_check_variable==boundary_condition_moved_wall) THEN
2452  CALL field_parameter_set_update_local_dof(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2453  & field_u_variable_type,field_values_set_type,local_ny, &
2454  & mesh_velocity_values(local_ny),err,error,*999)
2455  ELSE IF(boundary_condition_check_variable==boundary_condition_fixed_inlet) THEN
2456  CALL field_parameter_set_update_local_dof(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2457  & field_u_variable_type,field_values_set_type,local_ny, &
2458  & boundary_values(local_ny),err,error,*999)
2459  END IF
2460  ENDDO !deriv_idx
2461  ENDDO !node_idx
2462  ENDIF
2463  ENDIF
2464  ENDIF
2465  ENDDO !component_idx
2466  ENDIF
2467  ENDDO !variable_idx
2468  ELSE
2469  CALL flagerror("Boundary condition variable is not associated.",err,error,*999)
2470  END IF
2471  ELSE
2472  CALL flagerror("Boundary conditions are not associated.",err,error,*999)
2473  END IF
2474  ELSE
2475  CALL flagerror("Equations set is not associated.",err,error,*999)
2476  END IF
2477  ELSE
2478  CALL flagerror("Equations are not associated.",err,error,*999)
2479  END IF
2480  ELSE
2481  CALL flagerror("Solver equations are not associated.",err,error,*999)
2482  END IF
2483  CALL field_parameter_set_update_start(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2484  & field_values_set_type,err,error,*999)
2485  CALL field_parameter_set_update_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2486  & field_values_set_type,err,error,*999)
2487  END IF
2488  CASE(problem_ale_stokes_subtype)
2489  !Pre solve for the linear solver
2490  IF(solver%SOLVE_TYPE==solver_linear_type) THEN
2491  CALL write_string(general_output_type,"Mesh movement change boundary conditions... ",err,error,*999)
2492  solver_equations=>solver%SOLVER_EQUATIONS
2493  IF(ASSOCIATED(solver_equations)) THEN
2494  solver_mapping=>solver_equations%SOLVER_MAPPING
2495  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
2496  IF(ASSOCIATED(equations)) THEN
2497  equations_set=>equations%EQUATIONS_SET
2498  IF(ASSOCIATED(equations_set)) THEN
2499  boundary_conditions=>solver_equations%BOUNDARY_CONDITIONS
2500  IF(ASSOCIATED(boundary_conditions)) THEN
2501  CALL boundary_conditions_variable_get(boundary_conditions,equations_set%DEPENDENT%DEPENDENT_FIELD% &
2502  & variable_type_map(field_u_variable_type)%PTR,boundary_conditions_variable,err,error,*999)
2503  IF(ASSOCIATED(boundary_conditions_variable)) THEN
2504  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2505  & number_of_dimensions,err,error,*999)
2506  NULLIFY(boundary_values)
2507  CALL field_parameter_set_data_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
2508  & field_boundary_set_type,boundary_values,err,error,*999)
2509  CALL fluid_mechanics_io_read_boundary_conditions(solver_linear_type,boundary_values, &
2510  & number_of_dimensions,boundary_condition_moved_wall,control_loop%TIME_LOOP%INPUT_NUMBER, &
2511  & control_loop%TIME_LOOP%ITERATION_NUMBER,current_time,1.0_dp)
2512  DO variable_idx=1,equations_set%DEPENDENT%DEPENDENT_FIELD%NUMBER_OF_VARIABLES
2513  variable_type=equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
2514  field_variable=>equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
2515  IF(ASSOCIATED(field_variable)) THEN
2516  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
2517  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
2518  IF(ASSOCIATED(domain)) THEN
2519  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
2520  domain_nodes=>domain%TOPOLOGY%NODES
2521  IF(ASSOCIATED(domain_nodes)) THEN
2522  !Loop over the local nodes excluding the ghosts.
2523  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
2524  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
2525  !Default to version 1 of each node derivative
2526  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
2527  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
2528  boundary_condition_check_variable=boundary_conditions_variable% &
2529  & condition_types(local_ny)
2530  IF(boundary_condition_check_variable==boundary_condition_moved_wall) THEN
2531  CALL field_parameter_set_update_local_dof(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2532  & field_u_variable_type,field_values_set_type,local_ny, &
2533  & boundary_values(local_ny),err,error,*999)
2534  END IF
2535  END DO !deriv_idx
2536  ENDDO !node_idx
2537  ENDIF
2538  ENDIF
2539  ENDIF
2540  ENDDO !component_idx
2541  ENDIF
2542  ENDDO !variable_idx
2543  CALL field_parameter_set_data_restore(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
2544  & field_u_variable_type,field_boundary_set_type,boundary_values,err,error,*999)
2545 !\todo: This part should be read in out of a file eventually
2546  ELSE
2547  CALL flagerror("Boundary condition variable is not associated.",err,error,*999)
2548  END IF
2549  ELSE
2550  CALL flagerror("Boundary conditions are not associated.",err,error,*999)
2551  END IF
2552  ELSE
2553  CALL flagerror("Equations set is not associated.",err,error,*999)
2554  END IF
2555  ELSE
2556  CALL flagerror("Equations are not associated.",err,error,*999)
2557  END IF
2558  ELSE
2559  CALL flagerror("Solver equations are not associated.",err,error,*999)
2560  END IF
2561  CALL field_parameter_set_update_start(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2562  & field_values_set_type,err,error,*999)
2563  CALL field_parameter_set_update_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2564  & field_values_set_type,err,error,*999)
2565  !Pre solve for the dynamic solver
2566  ELSE IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
2567  CALL write_string(general_output_type,"Mesh movement change boundary conditions... ",err,error,*999)
2568  solver_equations=>solver%SOLVER_EQUATIONS
2569  IF(ASSOCIATED(solver_equations)) THEN
2570  solver_mapping=>solver_equations%SOLVER_MAPPING
2571  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
2572  IF(ASSOCIATED(equations)) THEN
2573  equations_set=>equations%EQUATIONS_SET
2574  IF(ASSOCIATED(equations_set)) THEN
2575  boundary_conditions=>solver_equations%BOUNDARY_CONDITIONS
2576  IF(ASSOCIATED(boundary_conditions)) THEN
2577  CALL boundary_conditions_variable_get(boundary_conditions,equations_set%DEPENDENT%DEPENDENT_FIELD% &
2578  & variable_type_map(field_u_variable_type)%PTR,boundary_conditions_variable,err,error,*999)
2579  IF(ASSOCIATED(boundary_conditions_variable)) THEN
2580  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2581  & number_of_dimensions,err,error,*999)
2582  NULLIFY(mesh_velocity_values)
2583  CALL field_parameter_set_data_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
2584  & field_mesh_velocity_set_type,mesh_velocity_values,err,error,*999)
2585  NULLIFY(boundary_values)
2586  CALL field_parameter_set_data_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
2587  & field_boundary_set_type,boundary_values,err,error,*999)
2588  CALL fluid_mechanics_io_read_boundary_conditions(solver_linear_type,boundary_values, &
2589  & number_of_dimensions,boundary_condition_fixed_inlet,control_loop%TIME_LOOP%INPUT_NUMBER, &
2590  & control_loop%TIME_LOOP%ITERATION_NUMBER,current_time,1.0_dp)
2591  DO variable_idx=1,equations_set%DEPENDENT%DEPENDENT_FIELD%NUMBER_OF_VARIABLES
2592  variable_type=equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
2593  field_variable=>equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
2594  IF(ASSOCIATED(field_variable)) THEN
2595  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
2596  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
2597  IF(ASSOCIATED(domain)) THEN
2598  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
2599  domain_nodes=>domain%TOPOLOGY%NODES
2600  IF(ASSOCIATED(domain_nodes)) THEN
2601  !Loop over the local nodes excluding the ghosts.
2602  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
2603  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
2604  !Default to version 1 of each node derivative
2605  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
2606  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
2607  displacement_value=0.0_dp
2608  boundary_condition_check_variable=boundary_conditions_variable% &
2609  & condition_types(local_ny)
2610  IF(boundary_condition_check_variable==boundary_condition_moved_wall) THEN
2611  CALL field_parameter_set_update_local_dof(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2612  & field_u_variable_type,field_values_set_type,local_ny, &
2613  & mesh_velocity_values(local_ny),err,error,*999)
2614  ELSE IF(boundary_condition_check_variable==boundary_condition_fixed_inlet) THEN
2615  CALL field_parameter_set_update_local_dof(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2616  & field_u_variable_type,field_values_set_type,local_ny, &
2617  & boundary_values(local_ny),err,error,*999)
2618  END IF
2619  END DO !deriv_idx
2620  ENDDO !node_idx
2621  ENDIF
2622  ENDIF
2623  ENDIF
2624  ENDDO !component_idx
2625  ENDIF
2626  ENDDO !variable_idx
2627  CALL field_parameter_set_data_restore(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
2628  & field_u_variable_type,field_mesh_velocity_set_type,mesh_velocity_values,err,error,*999)
2629  CALL field_parameter_set_data_restore(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
2630  & field_u_variable_type,field_boundary_set_type,boundary_values,err,error,*999)
2631  ELSE
2632  CALL flagerror("Boundary condition variable is not associated.",err,error,*999)
2633  END IF
2634  ELSE
2635  CALL flagerror("Boundary conditions are not associated.",err,error,*999)
2636  END IF
2637  ELSE
2638  CALL flagerror("Equations set is not associated.",err,error,*999)
2639  END IF
2640  ELSE
2641  CALL flagerror("Equations are not associated.",err,error,*999)
2642  END IF
2643  ELSE
2644  CALL flagerror("Solver equations are not associated.",err,error,*999)
2645  END IF
2646  CALL field_parameter_set_update_start(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2647  & field_values_set_type,err,error,*999)
2648  CALL field_parameter_set_update_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2649  & field_values_set_type,err,error,*999)
2650  END IF
2651  ! do nothing ???
2652  CASE DEFAULT
2653  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
2654  & " is not valid for a Stokes equation fluid type of a fluid mechanics problem class."
2655  CALL flagerror(local_error,err,error,*999)
2656  END SELECT
2657  ELSE
2658  CALL flagerror("Problem is not associated.",err,error,*999)
2659  ENDIF
2660  ELSE
2661  CALL flagerror("Solver is not associated.",err,error,*999)
2662  ENDIF
2663  ELSE
2664  CALL flagerror("Control loop is not associated.",err,error,*999)
2665  ENDIF
2666  exits("STOKES_PRE_SOLVE_UPDATE_BOUNDARY_CONDITIONS")
2667  RETURN
2668 999 errorsexits("STOKES_PRE_SOLVE_UPDATE_BOUNDARY_CONDITIONS",err,error)
2669  RETURN 1
2670  END SUBROUTINE stokes_pre_solve_update_boundary_conditions
2671 
2672  !
2673  !================================================================================================================================
2674  !
2676  SUBROUTINE stokes_pre_solve_ale_update_mesh(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
2677 
2678  !Argument variables
2679  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
2680  TYPE(solver_type), POINTER :: SOLVER
2681  INTEGER(INTG), INTENT(OUT) :: ERR
2682  TYPE(varying_string), INTENT(OUT) :: ERROR
2683  !Local Variables
2684  TYPE(solver_type), POINTER :: SOLVER_ALE_STOKES, SOLVER_LAPLACE
2685  TYPE(field_type), POINTER :: DEPENDENT_FIELD_LAPLACE, INDEPENDENT_FIELD_ALE_STOKES
2686  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS_LAPLACE, SOLVER_EQUATIONS_ALE_STOKES
2687  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING_LAPLACE, SOLVER_MAPPING_ALE_STOKES
2688  TYPE(equations_set_type), POINTER :: EQUATIONS_SET_LAPLACE, EQUATIONS_SET_ALE_STOKES
2689  TYPE(equations_type), POINTER :: EQUATIONS
2690  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
2691  TYPE(varying_string) :: LOCAL_ERROR
2692  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
2693  TYPE(domain_type), POINTER :: DOMAIN
2694  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
2695 
2696  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT,ALPHA
2697  REAL(DP), POINTER :: MESH_DISPLACEMENT_VALUES(:)
2698  INTEGER(INTG) :: I,NUMBER_OF_DIMENSIONS_LAPLACE,NUMBER_OF_DIMENSIONS_ALE_STOKES,GEOMETRIC_MESH_COMPONENT
2699  INTEGER(INTG) :: INPUT_TYPE,INPUT_OPTION,component_idx,deriv_idx,local_ny,node_idx,variable_idx,variable_type
2700 
2701  enters("STOKES_PRE_SOLVE_ALE_UPDATE_MESH",err,error,*999)
2702 
2703  IF(ASSOCIATED(control_loop)) THEN
2704  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
2705  NULLIFY(solver_laplace)
2706  NULLIFY(solver_ale_stokes)
2707  IF(ASSOCIATED(solver)) THEN
2708  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
2709  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
2710  CALL flagerror("Problem specification is not allocated.",err,error,*999)
2711  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
2712  CALL flagerror("Problem specification must have three entries for a Stokes problem.",err,error,*999)
2713  END IF
2714  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
2715  CASE(problem_static_stokes_subtype,problem_laplace_stokes_subtype)
2716  ! do nothing ???
2717  CASE(problem_transient_stokes_subtype)
2718  ! do nothing ???
2719  CASE(problem_pgm_stokes_subtype)
2720  !Update mesh within the dynamic solver
2721  IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
2722  !Get the independent field for the ALE Stokes problem
2723  CALL solvers_solver_get(solver%SOLVERS,1,solver_ale_stokes,err,error,*999)
2724  solver_equations_ale_stokes=>solver_ale_stokes%SOLVER_EQUATIONS
2725  IF(ASSOCIATED(solver_equations_ale_stokes)) THEN
2726  solver_mapping_ale_stokes=>solver_equations_ale_stokes%SOLVER_MAPPING
2727  IF(ASSOCIATED(solver_mapping_ale_stokes)) THEN
2728  equations_set_ale_stokes=>solver_mapping_ale_stokes%EQUATIONS_SETS(1)%PTR
2729  IF(ASSOCIATED(equations_set_ale_stokes)) THEN
2730  independent_field_ale_stokes=>equations_set_ale_stokes%INDEPENDENT%INDEPENDENT_FIELD
2731  ELSE
2732  CALL flagerror("ALE Stokes equations set is not associated.",err,error,*999)
2733  END IF
2734  !Get the data
2735  CALL field_number_of_components_get(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2736  & field_u_variable_type,number_of_dimensions_ale_stokes,err,error,*999)
2737 !\todo: Introduce flags set by the user (42/1 only for testings purpose)
2738  !Copy input to Stokes' independent field
2739  input_type=42
2740  input_option=1
2741  NULLIFY(mesh_displacement_values)
2742  CALL field_parameter_set_data_get(equations_set_ale_stokes%INDEPENDENT%INDEPENDENT_FIELD, &
2743  & field_u_variable_type,field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
2744  CALL fluid_mechanics_io_read_data(solver_linear_type,mesh_displacement_values, &
2745  & number_of_dimensions_ale_stokes,input_type,input_option,control_loop%TIME_LOOP%ITERATION_NUMBER,1.0_dp)
2746  CALL field_parameter_set_update_start(equations_set_ale_stokes%INDEPENDENT%INDEPENDENT_FIELD, &
2747  & field_u_variable_type,field_mesh_displacement_set_type,err,error,*999)
2748  CALL field_parameter_set_update_finish(equations_set_ale_stokes%INDEPENDENT%INDEPENDENT_FIELD, &
2749  & field_u_variable_type,field_mesh_displacement_set_type,err,error,*999)
2750  ELSE
2751  CALL flagerror("ALE Stokes solver mapping is not associated.",err,error,*999)
2752  END IF
2753  ELSE
2754  CALL flagerror("ALE Stokes solver equations are not associated.",err,error,*999)
2755  END IF
2756  !Use calculated values to update mesh
2757  CALL field_component_mesh_component_get(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2758  & field_u_variable_type,1,geometric_mesh_component,err,error,*999)
2759 ! CALL FIELD_PARAMETER_SET_DATA_GET(INDEPENDENT_FIELD_ALE_STOKES,FIELD_U_VARIABLE_TYPE, &
2760 ! & FIELD_MESH_DISPLACEMENT_SET_TYPE,MESH_DISPLACEMENT_VALUES,ERR,ERROR,*999)
2761  equations=>solver_mapping_ale_stokes%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
2762  IF(ASSOCIATED(equations)) THEN
2763  equations_mapping=>equations%EQUATIONS_MAPPING
2764  IF(ASSOCIATED(equations_mapping)) THEN
2765  DO variable_idx=1,equations_set_ale_stokes%DEPENDENT%DEPENDENT_FIELD%NUMBER_OF_VARIABLES
2766  variable_type=equations_set_ale_stokes%DEPENDENT%DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
2767  field_variable=>equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
2768  IF(ASSOCIATED(field_variable)) THEN
2769  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
2770  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
2771  IF(ASSOCIATED(domain)) THEN
2772  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
2773  domain_nodes=>domain%TOPOLOGY%NODES
2774  IF(ASSOCIATED(domain_nodes)) THEN
2775  !Loop over the local nodes excluding the ghosts.
2776  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
2777  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
2778  !Default to version 1 of each node derivative
2779  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
2780  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
2781  CALL field_parameter_set_add_local_dof(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2782  & field_u_variable_type,field_values_set_type,local_ny, &
2783  & mesh_displacement_values(local_ny),err,error,*999)
2784  ENDDO !deriv_idx
2785  ENDDO !node_idx
2786  ENDIF
2787  ENDIF
2788  ENDIF
2789  ENDDO !component_idx
2790  ENDIF
2791  ENDDO !variable_idx
2792  ELSE
2793  CALL flagerror("Equations mapping is not associated.",err,error,*999)
2794  END IF
2795  ELSE
2796  CALL flagerror("Equations are not associated.",err,error,*999)
2797  END IF
2798  CALL field_parameter_set_update_start(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2799  & field_u_variable_type,field_values_set_type,err,error,*999)
2800  CALL field_parameter_set_update_finish(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2801  & field_u_variable_type,field_values_set_type,err,error,*999)
2802  !Now use displacement values to calculate velocity values
2803  time_increment=control_loop%TIME_LOOP%TIME_INCREMENT
2804  alpha=1.0_dp/time_increment
2805  CALL field_parameter_sets_copy(independent_field_ale_stokes,field_u_variable_type, &
2806  & field_mesh_displacement_set_type,field_mesh_velocity_set_type,alpha,err,error,*999)
2807  ELSE
2808  CALL flagerror("Mesh motion calculation not successful for ALE problem.",err,error,*999)
2809  END IF
2810  CASE(problem_ale_stokes_subtype)
2811  !Update mesh within the dynamic solver
2812  IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
2813  IF(solver%DYNAMIC_SOLVER%ALE) THEN
2814  !Get the dependent field for the three component Laplace problem
2815  CALL solvers_solver_get(solver%SOLVERS,1,solver_laplace,err,error,*999)
2816  solver_equations_laplace=>solver_laplace%SOLVER_EQUATIONS
2817  IF(ASSOCIATED(solver_equations_laplace)) THEN
2818  solver_mapping_laplace=>solver_equations_laplace%SOLVER_MAPPING
2819  IF(ASSOCIATED(solver_mapping_laplace)) THEN
2820  equations_set_laplace=>solver_mapping_laplace%EQUATIONS_SETS(1)%PTR
2821  IF(ASSOCIATED(equations_set_laplace)) THEN
2822  dependent_field_laplace=>equations_set_laplace%DEPENDENT%DEPENDENT_FIELD
2823  ELSE
2824  CALL flagerror("Laplace equations set is not associated.",err,error,*999)
2825  END IF
2826  CALL field_number_of_components_get(equations_set_laplace%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2827  & number_of_dimensions_laplace,err,error,*999)
2828  ELSE
2829  CALL flagerror("Laplace solver mapping is not associated.",err,error,*999)
2830  END IF
2831  ELSE
2832  CALL flagerror("Laplace solver equations are not associated.",err,error,*999)
2833  END IF
2834  !Get the independent field for the ALE Stokes problem
2835  CALL solvers_solver_get(solver%SOLVERS,2,solver_ale_stokes,err,error,*999)
2836  solver_equations_ale_stokes=>solver_ale_stokes%SOLVER_EQUATIONS
2837  IF(ASSOCIATED(solver_equations_ale_stokes)) THEN
2838  solver_mapping_ale_stokes=>solver_equations_ale_stokes%SOLVER_MAPPING
2839  IF(ASSOCIATED(solver_mapping_ale_stokes)) THEN
2840  equations_set_ale_stokes=>solver_mapping_ale_stokes%EQUATIONS_SETS(1)%PTR
2841  IF(ASSOCIATED(equations_set_ale_stokes)) THEN
2842  independent_field_ale_stokes=>equations_set_ale_stokes%INDEPENDENT%INDEPENDENT_FIELD
2843  ELSE
2844  CALL flagerror("ALE Stokes equations set is not associated.",err,error,*999)
2845  END IF
2846  CALL field_number_of_components_get(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2847  & field_u_variable_type,number_of_dimensions_ale_stokes,err,error,*999)
2848  ELSE
2849  CALL flagerror("ALE Stokes solver mapping is not associated.",err,error,*999)
2850  END IF
2851  ELSE
2852  CALL flagerror("ALE Stokes solver equations are not associated.",err,error,*999)
2853  END IF
2854  !Copy result from Laplace mesh movement to Stokes' independent field
2855  IF(number_of_dimensions_ale_stokes==number_of_dimensions_laplace) THEN
2856  DO i=1,number_of_dimensions_ale_stokes
2857  CALL field_parameterstofieldparameterscopy(dependent_field_laplace, &
2858  & field_u_variable_type,field_values_set_type,i,independent_field_ale_stokes, &
2859  & field_u_variable_type,field_mesh_displacement_set_type,i,err,error,*999)
2860  END DO
2861  ELSE
2862  CALL flagerror("Dimension of Laplace and ALE Stokes equations set is not consistent.",err,error,*999)
2863  END IF
2864  !Use calculated values to update mesh
2865  CALL field_component_mesh_component_get(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2866  & field_u_variable_type,1,geometric_mesh_component,err,error,*999)
2867  NULLIFY(mesh_displacement_values)
2868  CALL field_parameter_set_data_get(independent_field_ale_stokes,field_u_variable_type, &
2869  & field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
2870  equations=>solver_mapping_laplace%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
2871  IF(ASSOCIATED(equations)) THEN
2872  equations_mapping=>equations%EQUATIONS_MAPPING
2873  IF(ASSOCIATED(equations_mapping)) THEN
2874  DO variable_idx=1,equations_set_ale_stokes%DEPENDENT%DEPENDENT_FIELD%NUMBER_OF_VARIABLES
2875  variable_type=equations_set_ale_stokes%DEPENDENT%DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
2876  field_variable=>equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
2877  IF(ASSOCIATED(field_variable)) THEN
2878  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
2879  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
2880  IF(ASSOCIATED(domain)) THEN
2881  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
2882  domain_nodes=>domain%TOPOLOGY%NODES
2883  IF(ASSOCIATED(domain_nodes)) THEN
2884  !Loop over the local nodes excluding the ghosts.
2885  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
2886  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
2887  !Default to version 1 of each node derivative
2888  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
2889  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
2890  CALL field_parameter_set_add_local_dof(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD, &
2891  & field_u_variable_type,field_values_set_type,local_ny, &
2892  & mesh_displacement_values(local_ny),err,error,*999)
2893  ENDDO !deriv_idx
2894  ENDDO !node_idx
2895  ENDIF
2896  ENDIF
2897  ENDIF
2898  ENDDO !component_idx
2899  ENDIF
2900  ENDDO !variable_idx
2901  ELSE
2902  CALL flagerror("Equations mapping is not associated.",err,error,*999)
2903  ENDIF
2904  CALL field_parameter_set_data_restore(independent_field_ale_stokes,field_u_variable_type, &
2905  & field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
2906  ELSE
2907  CALL flagerror("Equations are not associated.",err,error,*999)
2908  END IF
2909  CALL field_parameter_set_update_start(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2910  & field_values_set_type,err,error,*999)
2911  CALL field_parameter_set_update_finish(equations_set_ale_stokes%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2912  & field_values_set_type,err,error,*999)
2913  !Now use displacement values to calculate velocity values
2914  time_increment=control_loop%TIME_LOOP%TIME_INCREMENT
2915  alpha=1.0_dp/time_increment
2916  CALL field_parameter_sets_copy(independent_field_ale_stokes,field_u_variable_type, &
2917  & field_mesh_displacement_set_type,field_mesh_velocity_set_type,alpha,err,error,*999)
2918  ELSE
2919  CALL flagerror("Mesh motion calculation not successful for ALE problem.",err,error,*999)
2920  END IF
2921  ELSE
2922  CALL flagerror("Mesh update is not defined for non-dynamic problems.",err,error,*999)
2923  END IF
2924  CASE DEFAULT
2925  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
2926  & " is not valid for a Stokes equation fluid type of a fluid mechanics problem class."
2927  CALL flagerror(local_error,err,error,*999)
2928  END SELECT
2929  ELSE
2930  CALL flagerror("Problem is not associated.",err,error,*999)
2931  ENDIF
2932  ELSE
2933  CALL flagerror("Solver is not associated.",err,error,*999)
2934  ENDIF
2935  ELSE
2936  CALL flagerror("Control loop is not associated.",err,error,*999)
2937  ENDIF
2938  exits("STOKES_PRE_SOLVE_ALE_UPDATE_MESH")
2939  RETURN
2940 999 errorsexits("STOKES_PRE_SOLVE_ALE_UPDATE_MESH",err,error)
2941  RETURN 1
2942  END SUBROUTINE stokes_pre_solve_ale_update_mesh
2943 
2944  !
2945  !================================================================================================================================
2946  !
2948  SUBROUTINE stokes_pre_solve_ale_update_parameters(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
2949 
2950  !Argument variables
2951  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
2952  TYPE(solver_type), POINTER :: SOLVER
2953  INTEGER(INTG), INTENT(OUT) :: ERR
2954  TYPE(varying_string), INTENT(OUT) :: ERROR
2955  !Local Variables
2956  TYPE(field_type), POINTER :: INDEPENDENT_FIELD
2957  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
2958  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
2959  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2960  TYPE(equations_type), POINTER :: EQUATIONS
2961  TYPE(varying_string) :: LOCAL_ERROR
2962  TYPE(domain_type), POINTER :: DOMAIN
2963  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
2964  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
2965 
2966  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
2967  INTEGER(INTG) :: component_idx,node_idx,deriv_idx,local_ny,variable_idx,variable_type
2968  REAL(DP), POINTER :: MESH_STIFF_VALUES(:)
2969 
2970 
2971  enters("STOKES_PRE_SOLVE_ALE_UPDATE_PARAMETERS",err,error,*999)
2972 
2973  IF(ASSOCIATED(control_loop)) THEN
2974  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
2975  IF(ASSOCIATED(solver)) THEN
2976  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
2977  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
2978  CALL flagerror("Problem specification is not allocated.",err,error,*999)
2979  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
2980  CALL flagerror("Problem specification must have three entries for a Stokes problem.",err,error,*999)
2981  END IF
2982  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
2983  CASE(problem_static_stokes_subtype,problem_laplace_stokes_subtype)
2984  ! do nothing ???
2985  CASE(problem_transient_stokes_subtype)
2986  ! do nothing ???
2987  CASE(problem_ale_stokes_subtype)
2988  IF(solver%SOLVE_TYPE==solver_linear_type) THEN
2989  !Get the independent field for the ALE Stokes problem
2990  solver_equations=>solver%SOLVER_EQUATIONS
2991  IF(ASSOCIATED(solver_equations)) THEN
2992  solver_mapping=>solver_equations%SOLVER_MAPPING
2993  IF(ASSOCIATED(solver_mapping)) THEN
2994  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
2995  NULLIFY(mesh_stiff_values)
2996  CALL field_parameter_set_data_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
2997  & field_values_set_type,mesh_stiff_values,err,error,*999)
2998  IF(ASSOCIATED(equations_set)) THEN
2999  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
3000  IF(ASSOCIATED(equations)) THEN
3001  independent_field=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
3002  IF(ASSOCIATED(independent_field)) THEN
3003  DO variable_idx=1,equations_set%DEPENDENT%DEPENDENT_FIELD%NUMBER_OF_VARIABLES
3004  variable_type=equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
3005  field_variable=>equations_set%DEPENDENT%DEPENDENT_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
3006  IF(ASSOCIATED(field_variable)) THEN
3007  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
3008  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
3009  IF(ASSOCIATED(domain)) THEN
3010  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
3011  domain_nodes=>domain%TOPOLOGY%NODES
3012  IF(ASSOCIATED(domain_nodes)) THEN
3013  !Loop over the local nodes excluding the ghosts.
3014  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
3015  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
3016  !Default to version 1 of each node derivative
3017  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
3018  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
3019  ! !Calculation of K values dependent on current mesh topology
3020  mesh_stiff_values(local_ny)=1.0_dp
3021  CALL field_parameter_set_update_local_dof(equations_set%INDEPENDENT% &
3022  & independent_field,field_u_variable_type,field_values_set_type,local_ny, &
3023  & mesh_stiff_values(local_ny),err,error,*999)
3024  ENDDO !deriv_idx
3025  ENDDO !node_idx
3026  ENDIF
3027  ENDIF
3028  ENDIF
3029  ENDDO !component_idx
3030  ENDIF
3031  ENDDO !variable_idx
3032  ELSE
3033  CALL flagerror("Independent field is not associated.",err,error,*999)
3034  END IF
3035  ELSE
3036  CALL flagerror("Equations are not associated.",err,error,*999)
3037  END IF
3038  ELSE
3039  CALL flagerror("Equations set is not associated.",err,error,*999)
3040  ENDIF
3041  CALL field_parameter_set_data_restore(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
3042  & field_values_set_type,mesh_stiff_values,err,error,*999)
3043  ELSE
3044  CALL flagerror("Solver mapping is not associated.",err,error,*999)
3045  END IF
3046  ELSE
3047  CALL flagerror("Solver equations are not associated.",err,error,*999)
3048  END IF
3049  ELSE IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
3050  CALL flagerror("Mesh motion calculation not successful for ALE problem.",err,error,*999)
3051  END IF
3052  CASE DEFAULT
3053  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
3054  & " is not valid for a Stokes equation fluid type of a fluid mechanics problem class."
3055  CALL flagerror(local_error,err,error,*999)
3056  END SELECT
3057  ELSE
3058  CALL flagerror("Problem is not associated.",err,error,*999)
3059  ENDIF
3060  ELSE
3061  CALL flagerror("Solver is not associated.",err,error,*999)
3062  ENDIF
3063  ELSE
3064  CALL flagerror("Control loop is not associated.",err,error,*999)
3065  ENDIF
3066  exits("STOKES_PRE_SOLVE_ALE_UPDATE_PARAMETERS")
3067  RETURN
3068 999 errorsexits("STOKES_PRE_SOLVE_ALE_UPDATE_PARAMETERS",err,error)
3069  RETURN 1
3070  END SUBROUTINE stokes_pre_solve_ale_update_parameters
3071 
3072  !
3073  !================================================================================================================================
3074  !
3075 
3077  SUBROUTINE stokes_post_solve_output_data(CONTROL_LOOP,SOLVER,err,error,*)
3078 
3079  !Argument variables
3080  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
3081  TYPE(solver_type), POINTER :: SOLVER
3082  INTEGER(INTG), INTENT(OUT) :: ERR
3083  TYPE(varying_string), INTENT(OUT) :: ERROR
3084  !Local Variables
3085  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
3086  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
3087  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
3088  TYPE(fields_type), POINTER :: Fields
3089  TYPE(varying_string) :: localError,METHOD,FILENAME
3090 
3091  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
3092  INTEGER(INTG) :: EQUATIONS_SET_IDX,CURRENT_LOOP_ITERATION,OUTPUT_ITERATION_NUMBER,NUMBER_OF_DIMENSIONS
3093  LOGICAL :: EXPORT_FIELD
3094  CHARACTER(14) :: OUTPUT_FILE
3095 
3096  enters("STOKES_POST_SOLVE_OUTPUT_DATA",err,error,*999)
3097 
3098  NULLIFY(solver_equations)
3099  NULLIFY(solver_mapping)
3100  NULLIFY(equations_set)
3101 
3102  IF(ASSOCIATED(control_loop)) THEN
3103  IF(ASSOCIATED(solver)) THEN
3104  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
3105  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
3106  CALL flagerror("Problem specification is not allocated.",err,error,*999)
3107  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
3108  CALL flagerror("Problem specification must have three entries for a Stokes problem.",err,error,*999)
3109  END IF
3110  CALL system('mkdir -p ./output')
3111  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
3112  CASE(problem_static_stokes_subtype,problem_laplace_stokes_subtype)
3113  solver_equations=>solver%SOLVER_EQUATIONS
3114  IF(ASSOCIATED(solver_equations)) THEN
3115  solver_mapping=>solver_equations%SOLVER_MAPPING
3116  IF(ASSOCIATED(solver_mapping)) THEN
3117  !Make sure the equations sets are up to date
3118  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
3119  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%ptr
3120  filename="./output/"//"STATIC_SOLUTION"
3121  method="FORTRAN"
3122  IF(solver%output_Type>=solver_progress_output) THEN
3123  CALL write_string(general_output_type,"...",err,error,*999)
3124  CALL write_string(general_output_type,"Now export fields... ",err,error,*999)
3125  ENDIF
3126  fields=>equations_set%REGION%FIELDS
3127  CALL field_io_nodes_export(fields,filename,method,err,error,*999)
3128  CALL field_io_elements_export(fields,filename,method,err,error,*999)
3129  NULLIFY(fields)
3130  ENDDO
3131  ENDIF
3132  ENDIF
3133  CASE(problem_transient_stokes_subtype,problem_ale_stokes_subtype,problem_pgm_stokes_subtype)
3134  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
3135  solver_equations=>solver%SOLVER_EQUATIONS
3136  IF(ASSOCIATED(solver_equations)) THEN
3137  solver_mapping=>solver_equations%SOLVER_MAPPING
3138  IF(ASSOCIATED(solver_mapping)) THEN
3139  !Make sure the equations sets are up to date
3140  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
3141  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%ptr
3142  current_loop_iteration=control_loop%TIME_LOOP%ITERATION_NUMBER
3143  output_iteration_number=control_loop%TIME_LOOP%OUTPUT_NUMBER
3144  IF(output_iteration_number/=0) THEN
3145  IF(control_loop%TIME_LOOP%CURRENT_TIME<=control_loop%TIME_LOOP%STOP_TIME) THEN
3146  IF(current_loop_iteration<10) THEN
3147  WRITE(output_file,'("TIME_STEP_000",I0)') current_loop_iteration
3148  ELSE IF(current_loop_iteration<100) THEN
3149  WRITE(output_file,'("TIME_STEP_00",I0)') current_loop_iteration
3150  ELSE IF(current_loop_iteration<1000) THEN
3151  WRITE(output_file,'("TIME_STEP_0",I0)') current_loop_iteration
3152  ELSE IF(current_loop_iteration<10000) THEN
3153  WRITE(output_file,'("TIME_STEP_",I0)') current_loop_iteration
3154  END IF
3155  filename="./output/"//"MainTime_"//trim(numbertovstring(current_loop_iteration,"*",err,error))
3156  method="FORTRAN"
3157  IF(mod(current_loop_iteration,output_iteration_number)==0) THEN
3158  IF(control_loop%output_type >= control_loop_progress_output) THEN
3159  CALL write_string(general_output_type,"...",err,error,*999)
3160  CALL write_string(general_output_type,"Now export fields... ",err,error,*999)
3161  ENDIF
3162  fields=>equations_set%REGION%FIELDS
3163  CALL field_io_nodes_export(fields,filename,method,err,error,*999)
3164  CALL field_io_elements_export(fields,filename,method,err,error,*999)
3165  NULLIFY(fields)
3166  IF(control_loop%OUTPUT_TYPE >= control_loop_progress_output) THEN
3167  CALL write_string(general_output_type,filename,err,error,*999)
3168  CALL write_string(general_output_type,"...",err,error,*999)
3169  ENDIF
3170  END IF
3171  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
3172  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_navier_stokes_equation_two_dim_4.OR. &
3173  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_navier_stokes_equation_two_dim_5.OR. &
3174  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_navier_stokes_equation_three_dim_4.OR. &
3175  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_navier_stokes_equation_three_dim_5.OR. &
3176  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_navier_stokes_equation_three_dim_1) THEN
3177  CALL analyticanalysis_output(equations_set%DEPENDENT%DEPENDENT_FIELD,output_file,err,error,*999)
3178  ENDIF
3179  ENDIF
3180  ENDIF
3181  ENDIF
3182  ENDDO
3183  ENDIF
3184  ENDIF
3185  CASE DEFAULT
3186  localerror="Problem subtype "//trim(numbertovstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
3187  & " is not valid for a Stokes equation fluid type of a fluid mechanics problem class."
3188  CALL flagerror(localerror,err,error,*999)
3189  END SELECT
3190  ELSE
3191  CALL flagerror("Problem is not associated.",err,error,*999)
3192  ENDIF
3193  ELSE
3194  CALL flagerror("Solver is not associated.",err,error,*999)
3195  ENDIF
3196  ELSE
3197  CALL flagerror("Control loop is not associated.",err,error,*999)
3198  ENDIF
3199  exits("STOKES_POST_SOLVE_OUTPUT_DATA")
3200  RETURN
3201 999 errorsexits("STOKES_POST_SOLVE_OUTPUT_DATA",err,error)
3202  RETURN 1
3203  END SUBROUTINE stokes_post_solve_output_data
3204 
3205  !
3206  !================================================================================================================================
3207  !
3208 
3210  SUBROUTINE stokes_boundaryconditionsanalyticcalculate(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*)
3211 
3212  !Argument variables
3213  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
3214  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
3215  INTEGER(INTG), INTENT(OUT) :: ERR
3216  TYPE(varying_string), INTENT(OUT) :: ERROR
3217  !Local Variables
3218 !\todo: Reduce number of variables used
3219  INTEGER(INTG) :: component_idx,deriv_idx,dim_idx,local_ny,node_idx,NUMBER_OF_DIMENSIONS,variable_idx,variable_type,I,J,K
3220  INTEGER(INTG) :: number_of_nodes_xic(3),element_idx,en_idx,BOUND_COUNT,ANALYTIC_FUNCTION_TYPE,GLOBAL_DERIV_INDEX
3221  REAL(DP) :: VALUE,X(3),XI_COORDINATES(3)
3222 ! REAL(DP) :: BOUNDARY_TOLERANCE, BOUNDARY_X(3,2),MU_PARAM,L
3223  REAL(DP) :: T_COORDINATES(20,3),CURRENT_TIME,MU_PARAM,RHO_PARAM
3224  REAL(DP), POINTER :: GEOMETRIC_PARAMETERS(:)
3225  TYPE(domain_type), POINTER :: DOMAIN
3226  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
3227  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
3228  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
3229  TYPE(field_interpolated_point_ptr_type), POINTER :: INTERPOLATED_POINT(:)
3230  TYPE(field_interpolation_parameters_ptr_type), POINTER :: INTERPOLATION_PARAMETERS(:)
3231 ! TYPE(VARYING_STRING) :: LOCAL_ERROR
3232 
3233 ! ! ! !Temp variables
3234 ! ! ! INTEGER(INTG) :: number_of_element_nodes,temp_local_ny,temp_node_number,velocity_DOF_check,temp_local_node_number
3235 
3236  enters("Stokes_BoundaryConditionsAnalyticCalculate",err,error,*999)
3237 !\todo: Introduce user call to set parameters
3238  bound_count=0
3239 ! ! ! ! L=10.0_DP
3240  xi_coordinates(3)=0.0_dp
3241 ! BOUNDARY_TOLERANCE=0.000000001_DP
3242 ! ! ! BOUNDARY_X=0.0_DP
3243 ! ! ! T_COORDINATES=0.0_DP
3244 ! ! ! number_of_element_nodes=0
3245 ! ! ! temp_local_node_number=0
3246 ! ! ! temp_local_ny=0
3247 ! ! ! temp_node_number=0
3248 ! ! ! velocity_DOF_check=0
3249  IF(ASSOCIATED(equations_set)) THEN
3250  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
3251  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
3252  IF(ASSOCIATED(dependent_field)) THEN
3253  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
3254  IF(ASSOCIATED(geometric_field)) THEN
3255  NULLIFY(interpolation_parameters)
3256  NULLIFY(interpolated_point)
3257  CALL field_interpolation_parameters_initialise(geometric_field,interpolation_parameters,err,error,*999)
3258  CALL field_interpolated_points_initialise(interpolation_parameters,interpolated_point,err,error,*999)
3259  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
3260 ! ! ! !\todo: Check adjacent element calculation / use boundary node flag instead / didn't work for simplex
3261 ! ! ! IF(NUMBER_OF_DIMENSIONS==2) THEN
3262 ! ! ! BOUNDARY_X(1,1)=0.0_DP
3263 ! ! ! BOUNDARY_X(1,2)=10.0_DP
3264 ! ! ! BOUNDARY_X(2,1)=0.0_DP
3265 ! ! ! BOUNDARY_X(2,2)=10.0_DP
3266 ! ! ! ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
3267 ! ! ! BOUNDARY_X(1,1)=-5.0_DP
3268 ! ! ! BOUNDARY_X(1,2)=5.0_DP
3269 ! ! ! BOUNDARY_X(2,1)=-5.0_DP
3270 ! ! ! BOUNDARY_X(2,2)=5.0_DP
3271 ! ! ! BOUNDARY_X(3,1)=-5.0_DP
3272 ! ! ! BOUNDARY_X(3,2)=5.0_DP
3273 ! ! ! ENDIF
3274  NULLIFY(geometric_variable)
3275  CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
3276  NULLIFY(geometric_parameters)
3277  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,geometric_parameters, &
3278  & err,error,*999)
3279  IF(ASSOCIATED(boundary_conditions)) THEN
3280  DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES
3281  variable_type=dependent_field%VARIABLES(variable_idx)%VARIABLE_TYPE
3282  field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
3283  IF(ASSOCIATED(field_variable)) THEN
3284  CALL field_parameter_set_create(dependent_field,variable_type,field_analytic_values_set_type,err,error,*999)
3285  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
3286  bound_count=0
3287  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
3288  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
3289  IF(ASSOCIATED(domain)) THEN
3290  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
3291  domain_nodes=>domain%TOPOLOGY%NODES
3292  IF(ASSOCIATED(domain_nodes)) THEN
3293  !Loop over the local nodes excluding the ghosts.
3294  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
3295  element_idx=domain%topology%nodes%nodes(node_idx)%surrounding_elements(1)
3296  CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
3297  & interpolation_parameters(field_u_variable_type)%PTR,err,error,*999)
3298  en_idx=0
3299  xi_coordinates=0.0_dp
3300  number_of_nodes_xic(1)=domain%topology%elements%elements(element_idx)%basis%number_of_nodes_xic(1)
3301  number_of_nodes_xic(2)=domain%topology%elements%elements(element_idx)%basis%number_of_nodes_xic(2)
3302  IF(number_of_dimensions==3) THEN
3303  number_of_nodes_xic(3)=domain%topology%elements%elements(element_idx)%basis%number_of_nodes_xic(3)
3304  ELSE
3305  number_of_nodes_xic(3)=1
3306  ENDIF
3307  !\todo: Use boundary flag
3308  IF(domain%topology%elements%maximum_number_of_element_parameters==4.AND.number_of_dimensions==2 .OR. &
3309  & domain%topology%elements%maximum_number_of_element_parameters==9.OR. &
3310  & domain%topology%elements%maximum_number_of_element_parameters==16.OR. &
3311  & domain%topology%elements%maximum_number_of_element_parameters==8.OR. &
3312  & domain%topology%elements%maximum_number_of_element_parameters==27.OR. &
3313  & domain%topology%elements%maximum_number_of_element_parameters==64) THEN
3314  DO k=1,number_of_nodes_xic(3)
3315  DO j=1,number_of_nodes_xic(2)
3316  DO i=1,number_of_nodes_xic(1)
3317  en_idx=en_idx+1
3318  IF(domain%topology%elements%elements(element_idx)%element_nodes(en_idx)==node_idx) EXIT
3319  xi_coordinates(1)=xi_coordinates(1)+(1.0_dp/(number_of_nodes_xic(1)-1))
3320  ENDDO
3321  IF(domain%topology%elements%elements(element_idx)%element_nodes(en_idx)==node_idx) EXIT
3322  xi_coordinates(1)=0.0_dp
3323  xi_coordinates(2)=xi_coordinates(2)+(1.0_dp/(number_of_nodes_xic(2)-1))
3324  ENDDO
3325  IF(domain%topology%elements%elements(element_idx)%element_nodes(en_idx)==node_idx) EXIT
3326  xi_coordinates(1)=0.0_dp
3327  xi_coordinates(2)=0.0_dp
3328  IF(number_of_nodes_xic(3)/=1) THEN
3329  xi_coordinates(3)=xi_coordinates(3)+(1.0_dp/(number_of_nodes_xic(3)-1))
3330  ENDIF
3331  ENDDO
3332  CALL field_interpolate_xi(no_part_deriv,xi_coordinates, &
3333  & interpolated_point(field_u_variable_type)%PTR,err,error,*999)
3334  ELSE
3335  !\todo: Use boundary flag
3336  IF(domain%topology%elements%maximum_number_of_element_parameters==3) THEN
3337  t_coordinates(1,1:2)=[0.0_dp,1.0_dp]
3338  t_coordinates(2,1:2)=[1.0_dp,0.0_dp]
3339  t_coordinates(3,1:2)=[1.0_dp,1.0_dp]
3340  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==6) THEN
3341  t_coordinates(1,1:2)=[0.0_dp,1.0_dp]
3342  t_coordinates(2,1:2)=[1.0_dp,0.0_dp]
3343  t_coordinates(3,1:2)=[1.0_dp,1.0_dp]
3344  t_coordinates(4,1:2)=[0.5_dp,0.5_dp]
3345  t_coordinates(5,1:2)=[1.0_dp,0.5_dp]
3346  t_coordinates(6,1:2)=[0.5_dp,1.0_dp]
3347  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==10.AND. &
3348  & number_of_dimensions==2) THEN
3349  t_coordinates(1,1:2)=[0.0_dp,1.0_dp]
3350  t_coordinates(2,1:2)=[1.0_dp,0.0_dp]
3351  t_coordinates(3,1:2)=[1.0_dp,1.0_dp]
3352  t_coordinates(4,1:2)=[1.0_dp/3.0_dp,2.0_dp/3.0_dp]
3353  t_coordinates(5,1:2)=[2.0_dp/3.0_dp,1.0_dp/3.0_dp]
3354  t_coordinates(6,1:2)=[1.0_dp,1.0_dp/3.0_dp]
3355  t_coordinates(7,1:2)=[1.0_dp,2.0_dp/3.0_dp]
3356  t_coordinates(8,1:2)=[2.0_dp/3.0_dp,1.0_dp]
3357  t_coordinates(9,1:2)=[1.0_dp/3.0_dp,1.0_dp]
3358  t_coordinates(10,1:2)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp]
3359  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==4) THEN
3360  t_coordinates(1,1:3)=[0.0_dp,1.0_dp,1.0_dp]
3361  t_coordinates(2,1:3)=[1.0_dp,0.0_dp,1.0_dp]
3362  t_coordinates(3,1:3)=[1.0_dp,1.0_dp,0.0_dp]
3363  t_coordinates(4,1:3)=[1.0_dp,1.0_dp,1.0_dp]
3364  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==10.AND. &
3365  & number_of_dimensions==3) THEN
3366  t_coordinates(1,1:3)=[0.0_dp,1.0_dp,1.0_dp]
3367  t_coordinates(2,1:3)=[1.0_dp,0.0_dp,1.0_dp]
3368  t_coordinates(3,1:3)=[1.0_dp,1.0_dp,0.0_dp]
3369  t_coordinates(4,1:3)=[1.0_dp,1.0_dp,1.0_dp]
3370  t_coordinates(5,1:3)=[0.5_dp,0.5_dp,1.0_dp]
3371  t_coordinates(6,1:3)=[0.5_dp,1.0_dp,0.5_dp]
3372  t_coordinates(7,1:3)=[0.5_dp,1.0_dp,1.0_dp]
3373  t_coordinates(8,1:3)=[1.0_dp,0.5_dp,0.5_dp]
3374  t_coordinates(9,1:3)=[1.0_dp,1.0_dp,0.5_dp]
3375  t_coordinates(10,1:3)=[1.0_dp,0.5_dp,1.0_dp]
3376  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==20) THEN
3377  t_coordinates(1,1:3)=[0.0_dp,1.0_dp,1.0_dp]
3378  t_coordinates(2,1:3)=[1.0_dp,0.0_dp,1.0_dp]
3379  t_coordinates(3,1:3)=[1.0_dp,1.0_dp,0.0_dp]
3380  t_coordinates(4,1:3)=[1.0_dp,1.0_dp,1.0_dp]
3381  t_coordinates(5,1:3)=[1.0_dp/3.0_dp,2.0_dp/3.0_dp,1.0_dp]
3382  t_coordinates(6,1:3)=[2.0_dp/3.0_dp,1.0_dp/3.0_dp,1.0_dp]
3383  t_coordinates(7,1:3)=[1.0_dp/3.0_dp,1.0_dp,2.0_dp/3.0_dp]
3384  t_coordinates(8,1:3)=[2.0_dp/3.0_dp,1.0_dp,1.0_dp/3.0_dp]
3385  t_coordinates(9,1:3)=[1.0_dp/3.0_dp,1.0_dp,1.0_dp]
3386  t_coordinates(10,1:3)=[2.0_dp/3.0_dp,1.0_dp,1.0_dp]
3387  t_coordinates(11,1:3)=[1.0_dp,1.0_dp/3.0_dp,2.0_dp/3.0_dp]
3388  t_coordinates(12,1:3)=[1.0_dp,2.0_dp/3.0_dp,1.0_dp/3.0_dp]
3389  t_coordinates(13,1:3)=[1.0_dp,1.0_dp,1.0_dp/3.0_dp]
3390  t_coordinates(14,1:3)=[1.0_dp,1.0_dp,2.0_dp/3.0_dp]
3391  t_coordinates(15,1:3)=[1.0_dp,1.0_dp/3.0_dp,1.0_dp]
3392  t_coordinates(16,1:3)=[1.0_dp,2.0_dp/3.0_dp,1.0_dp]
3393  t_coordinates(17,1:3)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp,2.0_dp/3.0_dp]
3394  t_coordinates(18,1:3)=[2.0_dp/3.0_dp,2.0_dp/3.0_dp,1.0_dp]
3395  t_coordinates(19,1:3)=[2.0_dp/3.0_dp,1.0_dp,2.0_dp/3.0_dp]
3396  t_coordinates(20,1:3)=[1.0_dp,2.0_dp/3.0_dp,2.0_dp/3.0_dp]
3397  ENDIF
3398  DO k=1,domain%topology%elements%maximum_number_of_element_parameters
3399  IF(domain%topology%elements%elements(element_idx)%element_nodes(k)==node_idx) EXIT
3400  ENDDO
3401  IF(number_of_dimensions==2) THEN
3402  CALL field_interpolate_xi(no_part_deriv,t_coordinates(k,1:2), &
3403  & interpolated_point(field_u_variable_type)%PTR,err,error,*999)
3404  ELSE IF(number_of_dimensions==3) THEN
3405  CALL field_interpolate_xi(no_part_deriv,t_coordinates(k,1:3), &
3406  & interpolated_point(field_u_variable_type)%PTR,err,error,*999)
3407  ENDIF
3408  ENDIF
3409  x=0.0_dp
3410  DO dim_idx=1,number_of_dimensions
3411  x(dim_idx)=interpolated_point(field_u_variable_type)%PTR%VALUES(dim_idx,1)
3412  ENDDO !dim_idx
3413 
3414  !Loop over the derivatives
3415  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
3416  analytic_function_type=equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE
3417  global_deriv_index=domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX
3418  current_time=0.0_dp
3419  materials_field=>equations_set%MATERIALS%MATERIALS_FIELD
3420  !Define MU_PARAM, density=1
3421  mu_param=materials_field%variables(1)%parameter_sets%parameter_sets(1)%ptr% &
3422  & parameters%cmiss%data_dp(1)
3423  !Define RHO_PARAM, density=2
3424  IF(analytic_function_type==equations_set_stokes_equation_two_dim_4.OR. &
3425  & analytic_function_type==equations_set_stokes_equation_two_dim_5.OR. &
3426  & analytic_function_type==equations_set_stokes_equation_three_dim_4.OR. &
3427  & analytic_function_type==equations_set_stokes_equation_three_dim_5) THEN
3428  rho_param=materials_field%variables(1)%parameter_sets%parameter_sets(1)%ptr% &
3429  & parameters%cmiss%data_dp(2)
3430  ELSE
3431  rho_param=0.0_dp
3432  ENDIF
3433  CALL stokes_equation_analytic_functions(VALUE,x,mu_param,rho_param,current_time,variable_type, &
3434  & global_deriv_index,analytic_function_type,number_of_dimensions, &
3435  & field_variable%NUMBER_OF_COMPONENTS,component_idx,err,error,*999)
3436  !Default to version 1 of each node derivative
3437  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
3438  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
3439  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
3440  & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
3441  IF(variable_type==field_u_variable_type) THEN
3442  ! \todo: This part should work even for simplex elements as soon as adjacent element calculation has been fixed
3443  IF(domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN
3444  !If we are a boundary node then set the analytic value on the boundary
3445  IF(component_idx<=number_of_dimensions) THEN
3446  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field,variable_type, &
3447  & local_ny,boundary_condition_fixed,VALUE,err,error,*999)
3448  bound_count=bound_count+1
3449  ELSE
3450  ! \todo: This is just a workaround for linear pressure fields in simplex element components
3451  IF(domain%topology%elements%maximum_number_of_element_parameters==3) THEN
3452  IF(analytic_function_type==equations_set_stokes_equation_two_dim_1.OR. &
3453  & analytic_function_type==equations_set_stokes_equation_two_dim_2.OR. &
3454  & analytic_function_type==equations_set_stokes_equation_two_dim_3.OR. &
3455  & analytic_function_type==equations_set_stokes_equation_two_dim_4.OR. &
3456  & analytic_function_type==equations_set_stokes_equation_two_dim_5) THEN
3457  IF(-0.001_dp<x(1).AND.x(1)<0.001_dp.AND.-0.001_dp<x(2).AND.x(2)<0.001_dp.OR. &
3458  & 10.0_dp-0.001_dp<x(1).AND.x(1)<10.0_dp+0.001_dp.AND.-0.001_dp<x(2).AND. &
3459  & x(2)<0.001_dp.OR. &
3460  & 10.0_dp-0.001_dp<x(1).AND.x(1)<10.0_dp+0.001_dp.AND.10.0_dp-0.001_dp<x(2).AND. &
3461  & x(2)<10.0_dp+0.001_dp.OR. &
3462  & -0.001_dp<x(1).AND.x(1)<0.001_dp.AND.10.0_dp-0.001_dp<x(2).AND. &
3463  & x(2)<10.0_dp+0.001_dp) THEN
3464  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field, &
3465  & variable_type,local_ny,boundary_condition_fixed,VALUE,err,error,*999)
3466  bound_count=bound_count+1
3467  ENDIF
3468  ENDIF
3469  ELSE IF(domain%topology%elements%maximum_number_of_element_parameters==4.AND. &
3470  & number_of_dimensions==3) THEN
3471  IF(analytic_function_type==equations_set_stokes_equation_three_dim_1.OR. &
3472  & analytic_function_type==equations_set_stokes_equation_three_dim_2.OR. &
3473  & analytic_function_type==equations_set_stokes_equation_three_dim_3.OR. &
3474  & analytic_function_type==equations_set_stokes_equation_three_dim_4.OR. &
3475  & analytic_function_type==equations_set_stokes_equation_three_dim_5) THEN
3476  IF(-5.0_dp-0.001_dp<x(1).AND.x(1)<-5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(2).AND. &
3477  & x(2)<-5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(3).AND.x(3)<-5.0_dp+0.001_dp.OR. &
3478  & -5.0_dp-0.001_dp<x(1).AND.x(1)<-5.0_dp+0.001_dp.AND.5.0_dp-0.001_dp<x(2).AND. &
3479  & x(2)<5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(3).AND.x(3)<-5.0_dp+0.001_dp.OR. &
3480  & 5.0_dp-0.001_dp<x(1).AND.x(1)<5.0_dp+0.001_dp.AND.5.0_dp-0.001_dp<x(2).AND. &
3481  & x(2)<5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(3).AND.x(3)<-5.0_dp+0.001_dp.OR. &
3482  & 5.0_dp-0.001_dp<x(1).AND.x(1)<5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(2).AND. &
3483  & x(2)<-5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(3).AND.x(3)<-5.0_dp+0.001_dp.OR. &
3484  & -5.0_dp-0.001_dp<x(1).AND.x(1)<-5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(2).AND. &
3485  & x(2)<-5.0_dp+0.001_dp.AND.5.0_dp-0.001_dp<x(3).AND.x(3)<5.0_dp+0.001_dp.OR. &
3486  & -5.0_dp-0.001_dp<x(1).AND.x(1)<-5.0_dp+0.001_dp.AND.5.0_dp-0.001_dp<x(2).AND. &
3487  & x(2)<5.0_dp+0.001_dp.AND.5.0_dp-0.001_dp<x(3).AND.x(3)<5.0_dp+0.001_dp.OR. &
3488  & 5.0_dp-0.001_dp<x(1).AND.x(1)<5.0_dp+0.001_dp.AND.5.0_dp-0.001_dp<x(2).AND. &
3489  & x(2)<5.0_dp+0.001_dp.AND.5.0_dp-0.001_dp<x(3).AND.x(3)<5.0_dp+0.001_dp.OR. &
3490  & 5.0_dp-0.001_dp<x(1).AND.x(1)<5.0_dp+0.001_dp.AND.-5.0_dp-0.001_dp<x(2).AND. &
3491  & x(2)<-5.0_dp+ 0.001_dp.AND.5.0_dp-0.001_dp<x(3).AND.x(3)<5.0_dp+0.001_dp) THEN
3492  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field, &
3493  & variable_type,local_ny,boundary_condition_fixed,VALUE,err,error,*999)
3494  bound_count=bound_count+1
3495  ENDIF
3496  ENDIF
3497  ! \todo: This is how it should be if adjacent elements would be working
3498  ELSE IF(bound_count==0) THEN
3499  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field,variable_type, &
3500  & local_ny,boundary_condition_fixed,VALUE,err,error,*999)
3501  bound_count=bound_count+1
3502  ENDIF
3503 
3504 
3505  ENDIF
3506  ELSE
3507  IF(component_idx<=number_of_dimensions) THEN
3508  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
3509  & field_values_set_type,local_ny,VALUE,err,error,*999)
3510  ENDIF
3511  ENDIF
3512  ! \todo: Use boundary node flag
3513  ! ! ! !If we are a boundary node then set the analytic value on the boundary
3514  ! ! ! IF(NUMBER_OF_DIMENSIONS==2) THEN
3515  ! ! ! IF(X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.OR. &
3516  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.OR. &
3517  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.OR. &
3518  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE) THEN
3519  ! ! ! IF(component_idx<=NUMBER_OF_DIMENSIONS) THEN
3520  ! ! ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3521  ! ! ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3522  ! ! ! BOUND_COUNT=BOUND_COUNT+1
3523  ! ! ! !Apply boundary conditions check for pressure nodes
3524  ! ! ! ELSE IF(component_idx>NUMBER_OF_DIMENSIONS) THEN
3525  ! ! ! IF(DOMAIN%topology%elements%maximum_number_of_element_parameters==4) THEN
3526  ! ! ! IF(X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.AND. &
3527  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE) &
3528  ! ! ! & THEN
3529  ! ! ! ! Commented out for testing purposes
3530  ! ! ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3531  ! ! ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3532  ! ! ! BOUND_COUNT=BOUND_COUNT+1
3533  ! ! ! ENDIF
3534  ! ! ! ENDIF
3535  ! ! ! !\todo: Again, ...
3536  ! ! ! IF(DOMAIN%topology%elements%maximum_number_of_element_parameters==3.OR. &
3537  ! ! ! & DOMAIN%topology%elements%maximum_number_of_element_parameters==6.OR. &
3538  ! ! ! & DOMAIN%topology%elements%maximum_number_of_element_parameters==10) THEN
3539  ! ! ! IF(X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.AND. &
3540  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.OR. &
3541  ! ! ! & X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.AND.&
3542  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE.OR. &
3543  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.AND.&
3544  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.OR. &
3545  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.AND.&
3546  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE) &
3547  ! ! ! & THEN
3548  ! ! ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3549  ! ! ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3550  ! ! ! BOUND_COUNT=BOUND_COUNT+1
3551  ! ! ! ENDIF
3552  ! ! ! ENDIF
3553  ! ! ! ENDIF
3554  ! ! ! ENDIF
3555  ! ! ! IF(component_idx<=NUMBER_OF_DIMENSIONS+1) THEN
3556  ! ! ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_DOF(DEPENDENT_FIELD,variable_type, &
3557  ! ! ! & FIELD_VALUES_SET_TYPE,local_ny,VALUE,ERR,ERROR,*999)
3558  ! ! ! ENDIF
3559  ! ! ! ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
3560  ! ! ! IF(X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.OR. &
3561  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.OR. &
3562  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.OR. &
3563  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE.OR. &
3564  ! ! ! & X(3)<BOUNDARY_X(3,1)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,1)-BOUNDARY_TOLERANCE.OR. &
3565  ! ! ! & X(3)<BOUNDARY_X(3,2)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,2)-BOUNDARY_TOLERANCE) THEN
3566  ! ! ! IF(component_idx<=NUMBER_OF_DIMENSIONS) THEN
3567  ! ! ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3568  ! ! ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3569  ! ! ! BOUND_COUNT=BOUND_COUNT+1
3570  ! ! ! !Apply boundary conditions check for pressure nodes
3571  ! ! ! ELSE IF(component_idx>NUMBER_OF_DIMENSIONS) THEN
3572  ! ! ! IF(DOMAIN%topology%elements%maximum_number_of_element_parameters==4.OR. &
3573  ! ! ! & DOMAIN%topology%elements%maximum_number_of_element_parameters==10.OR. &
3574  ! ! ! & DOMAIN%topology%elements%maximum_number_of_element_parameters==20) THEN
3575  ! ! ! IF(X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.AND. &
3576  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.AND. &
3577  ! ! ! & X(3)<BOUNDARY_X(3,1)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,1)-BOUNDARY_TOLERANCE.OR. &
3578  ! ! ! & X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.AND. &
3579  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.AND. &
3580  ! ! ! & X(3)<BOUNDARY_X(3,2)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,2)-BOUNDARY_TOLERANCE.OR. &
3581  ! ! ! & X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.AND. &
3582  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE.AND. &
3583  ! ! ! & X(3)<BOUNDARY_X(3,1)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,1)-BOUNDARY_TOLERANCE.OR. &
3584  ! ! ! & X(1)<BOUNDARY_X(1,1)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,1)-BOUNDARY_TOLERANCE.AND. &
3585  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE.AND. &
3586  ! ! ! & X(3)<BOUNDARY_X(3,2)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,2)-BOUNDARY_TOLERANCE.OR. &
3587  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.AND. &
3588  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.AND. &
3589  ! ! ! & X(3)<BOUNDARY_X(3,1)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,1)-BOUNDARY_TOLERANCE.OR. &
3590  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.AND. &
3591  ! ! ! & X(2)<BOUNDARY_X(2,1)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,1)-BOUNDARY_TOLERANCE.AND. &
3592  ! ! ! & X(3)<BOUNDARY_X(3,2)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,2)-BOUNDARY_TOLERANCE.OR. &
3593  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.AND. &
3594  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE.AND. &
3595  ! ! ! & X(3)<BOUNDARY_X(3,1)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,1)-BOUNDARY_TOLERANCE.OR. &
3596  ! ! ! & X(1)<BOUNDARY_X(1,2)+BOUNDARY_TOLERANCE.AND.X(1)>BOUNDARY_X(1,2)-BOUNDARY_TOLERANCE.AND. &
3597  ! ! ! & X(2)<BOUNDARY_X(2,2)+BOUNDARY_TOLERANCE.AND.X(2)>BOUNDARY_X(2,2)-BOUNDARY_TOLERANCE.AND. &
3598  ! ! ! & X(3)<BOUNDARY_X(3,2)+BOUNDARY_TOLERANCE.AND.X(3)>BOUNDARY_X(3,2)-BOUNDARY_TOLERANCE) THEN
3599  ! ! ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3600  ! ! ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3601  ! ! ! BOUND_COUNT=BOUND_COUNT+1
3602  ! ! ! ENDIF
3603  ! ! ! ENDIF
3604  ! ! ! ENDIF
3605  ! ! ! ELSE
3606  ! ! ! IF(component_idx<=NUMBER_OF_DIMENSIONS+1) THEN
3607  ! ! ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_DOF(DEPENDENT_FIELD,variable_type, &
3608  ! ! ! & FIELD_VALUES_SET_TYPE,local_ny,VALUE,ERR,ERROR,*999)
3609  ! ! ! ENDIF
3610  ! ! ! ENDIF
3611  ! ! ! ENDIF
3612  ENDIF
3613  ENDDO !deriv_idx
3614  ENDDO !node_idx
3615  ELSE
3616  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
3617  ENDIF
3618  ELSE
3619  CALL flagerror("Domain topology is not associated.",err,error,*999)
3620  ENDIF
3621  ELSE
3622  CALL flagerror("Domain is not associated.",err,error,*999)
3623  ENDIF
3624  ELSE
3625  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
3626  ENDIF
3627  WRITE(*,*)'NUMBER OF BOUNDARIES SET ',bound_count
3628  ENDDO !component_idx
3629  CALL field_parameter_set_update_start(dependent_field,variable_type,field_analytic_values_set_type, &
3630  & err,error,*999)
3631  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_analytic_values_set_type, &
3632  & err,error,*999)
3633  CALL field_parameter_set_update_start(dependent_field,variable_type,field_values_set_type, &
3634  & err,error,*999)
3635  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_values_set_type, &
3636  & err,error,*999)
3637  ELSE
3638  CALL flagerror("Field variable is not associated.",err,error,*999)
3639  ENDIF
3640  ENDDO !variable_idx
3641  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
3642  & geometric_parameters,err,error,*999)
3643  CALL field_interpolated_points_finalise(interpolated_point,err,error,*999)
3644  CALL field_interpolation_parameters_finalise(interpolation_parameters,err,error,*999)
3645  ELSE
3646  CALL flagerror("Boundary conditions is not associated.",err,error,*999)
3647  ENDIF
3648  ELSE
3649  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
3650  ENDIF
3651  ELSE
3652  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
3653  ENDIF
3654  ELSE
3655  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
3656  ENDIF
3657  ELSE
3658  CALL flagerror("Equations set is not associated.",err,error,*999)
3659  ENDIF
3660 
3661  exits("Stokes_BoundaryConditionsAnalyticCalculate")
3662  RETURN
3663 999 errorsexits("Stokes_BoundaryConditionsAnalyticCalculate",err,error)
3664  RETURN 1
3665  END SUBROUTINE stokes_boundaryconditionsanalyticcalculate
3666 
3667  !
3668  !================================================================================================================================
3669  !
3671  SUBROUTINE stokes_equation_analytic_functions(VALUE,X,MU_PARAM,RHO_PARAM,CURRENT_TIME,VARIABLE_TYPE, &
3672  & global_deriv_index,analytic_function_type,number_of_dimensions,number_of_components,component_idx,err,error,*)
3673 
3674  !Argument variables
3675  INTEGER(INTG), INTENT(OUT) :: ERR
3676  TYPE(varying_string), INTENT(OUT) :: ERROR
3677  REAL(DP), INTENT(OUT) :: VALUE
3678  REAL(DP) :: MU_PARAM,RHO_PARAM
3679  REAL(DP), INTENT(IN) :: CURRENT_TIME
3680  REAL(DP), INTENT(IN), DIMENSION(3) :: X
3681  INTEGER(INTG), INTENT(IN) :: NUMBER_OF_DIMENSIONS,NUMBER_OF_COMPONENTS,COMPONENT_IDX
3682  !Local variables
3683  TYPE(varying_string) :: LOCAL_ERROR
3684  INTEGER(INTG) :: variable_type,GLOBAL_DERIV_INDEX,ANALYTIC_FUNCTION_TYPE
3685  !TYPE(DOMAIN_TYPE), POINTER :: DOMAIN
3686  !TYPE(DOMAIN_NODES_TYPE), POINTER :: DOMAIN_NODES
3687  REAL(DP) :: INTERNAL_TIME
3688 
3689  enters("STOKES_EQUATION_ANALYTIC_FUNCTIONS",err,error,*999)
3690 
3691 !\todo: Introduce user-defined or default values instead for density and viscosity
3692  internal_time=current_time
3693  SELECT CASE(analytic_function_type)
3694  CASE(equations_set_stokes_equation_two_dim_1)
3695  IF(number_of_dimensions==2.AND.number_of_components==3) THEN
3696  !Polynomial function
3697  SELECT CASE(variable_type)
3698  CASE(field_u_variable_type)
3699  SELECT CASE(global_deriv_index)
3700  CASE(no_global_deriv)
3701  IF(component_idx==1) THEN
3702  !calculate u
3703  VALUE=x(2)**2/10.0_dp**2
3704  ELSE IF(component_idx==2) THEN
3705  !calculate v
3706  VALUE=x(1)**2/10.0_dp**2
3707  ELSE IF(component_idx==3) THEN
3708  !calculate p
3709  VALUE=2.0_dp*mu_param/10.0_dp**2*x(1)
3710  ELSE
3711  CALL flagerror("Not implemented.",err,error,*999)
3712  ENDIF
3713  CASE(global_deriv_s1)
3714  CALL flagerror("Not implemented.",err,error,*999)
3715  CASE(global_deriv_s2)
3716  CALL flagerror("Not implemented.",err,error,*999)
3717  CASE(global_deriv_s1_s2)
3718  CALL flagerror("Not implemented.",err,error,*999)
3719  CASE DEFAULT
3720  local_error="The global derivative index of "//trim(number_to_vstring( &
3721  & global_deriv_index,"*",err,error))// &
3722  & " is invalid."
3723  CALL flagerror(local_error,err,error,*999)
3724  END SELECT
3725  CASE(field_deludeln_variable_type)
3726  SELECT CASE(global_deriv_index)
3727  CASE(no_global_deriv)
3728  VALUE= 0.0_dp
3729  CASE(global_deriv_s1)
3730  CALL flagerror("Not implemented.",err,error,*999)
3731  CASE(global_deriv_s2)
3732  CALL flagerror("Not implemented.",err,error,*999)
3733  CASE(global_deriv_s1_s2)
3734  CALL flagerror("Not implemented.",err,error,*999)
3735  CASE DEFAULT
3736  local_error="The global derivative index of "//trim(number_to_vstring( &
3737  & global_deriv_index,"*",err,error))// &
3738  & " is invalid."
3739  CALL flagerror(local_error,err,error,*999)
3740  END SELECT
3741  CASE DEFAULT
3742  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
3743  & " is invalid."
3744  CALL flagerror(local_error,err,error,*999)
3745  END SELECT
3746  ELSE
3747  local_error="The number of components does not correspond to the number of dimensions."
3748  CALL flagerror(local_error,err,error,*999)
3749  ENDIF
3750  CASE(equations_set_stokes_equation_two_dim_2)
3751  IF(number_of_dimensions==2.AND.number_of_components==3) THEN
3752  !Exponential function
3753  SELECT CASE(variable_type)
3754  CASE(field_u_variable_type)
3755  SELECT CASE(global_deriv_index)
3756  CASE(no_global_deriv)
3757  IF(component_idx==1) THEN
3758  !calculate u
3759  VALUE= exp((x(1)-x(2))/10.0_dp)
3760  ELSE IF(component_idx==2) THEN
3761  !calculate v
3762  VALUE= exp((x(1)-x(2))/10.0_dp)
3763  ELSE IF(component_idx==3) THEN
3764  !calculate p
3765  VALUE= 2.0_dp*mu_param/10.0_dp*exp((x(1)-x(2))/10.0_dp)
3766  ELSE
3767  CALL flagerror("Not implemented.",err,error,*999)
3768  ENDIF
3769  CASE(global_deriv_s1)
3770  CALL flagerror("Not implemented.",err,error,*999)
3771  CASE(global_deriv_s2)
3772  CALL flagerror("Not implemented.",err,error,*999)
3773  CASE(global_deriv_s1_s2)
3774  CALL flagerror("Not implemented.",err,error,*999)
3775  CASE DEFAULT
3776  local_error="The global derivative index of "//trim(number_to_vstring( &
3777  & global_deriv_index,"*",err,error))// &
3778  & " is invalid."
3779  CALL flagerror(local_error,err,error,*999)
3780  END SELECT
3781  CASE(field_deludeln_variable_type)
3782  SELECT CASE(global_deriv_index)
3783  CASE(no_global_deriv)
3784  IF(component_idx==1) THEN
3785  !calculate u
3786  VALUE= 0.0_dp
3787  ELSE IF(component_idx==2) THEN
3788  !calculate v
3789  VALUE= 0.0_dp
3790  ELSE IF(component_idx==3) THEN
3791  !calculate p
3792  VALUE= 0.0_dp
3793  ELSE
3794  CALL flagerror("Not implemented.",err,error,*999)
3795  ENDIF
3796  CASE(global_deriv_s1)
3797  CALL flagerror("Not implemented.",err,error,*999)
3798  CASE(global_deriv_s2)
3799  CALL flagerror("Not implemented.",err,error,*999)
3800  CASE(global_deriv_s1_s2)
3801  CALL flagerror("Not implemented.",err,error,*999)
3802  CASE DEFAULT
3803  local_error="The global derivative index of "//trim(number_to_vstring( &
3804  & global_deriv_index,"*",err,error))// &
3805  & " is invalid."
3806  CALL flagerror(local_error,err,error,*999)
3807  END SELECT
3808  CASE DEFAULT
3809  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
3810  & " is invalid."
3811  CALL flagerror(local_error,err,error,*999)
3812  END SELECT
3813  ELSE
3814  local_error="The number of components does not correspond to the number of dimensions."
3815  CALL flagerror(local_error,err,error,*999)
3816  ENDIF
3817  CASE(equations_set_stokes_equation_two_dim_3)
3818  IF(number_of_dimensions==2.AND.number_of_components==3) THEN
3819  !Sine and cosine functions
3820  SELECT CASE(variable_type)
3821  CASE(field_u_variable_type)
3822  SELECT CASE(global_deriv_index)
3823  CASE(no_global_deriv)
3824  IF(component_idx==1) THEN
3825  !calculate u
3826  VALUE=sin(2.0_dp*pi*x(1)/10.0_dp)*sin(2.0_dp*pi*x(2)/10.0_dp)
3827  ELSE IF(component_idx==2) THEN
3828  !calculate v
3829  VALUE=cos(2.0_dp*pi*x(1)/10.0_dp)*cos(2.0_dp*pi*x(2)/10.0_dp)
3830  ELSE IF(component_idx==3) THEN
3831  !calculate p
3832  VALUE=4.0_dp*mu_param*pi/10.0_dp*sin(2.0_dp*pi*x(2)/10.0_dp)*cos(2.0_dp*pi*x(1)/10.0_dp)
3833  ELSE
3834  CALL flagerror("Not implemented.",err,error,*999)
3835  ENDIF
3836  CASE(global_deriv_s1)
3837  CALL flagerror("Not implemented.",err,error,*999)
3838  CASE(global_deriv_s2)
3839  CALL flagerror("Not implemented.",err,error,*999)
3840  CASE(global_deriv_s1_s2)
3841  CALL flagerror("Not implemented.",err,error,*999)
3842  CASE DEFAULT
3843  local_error="The global derivative index of "//trim(number_to_vstring( &
3844  & global_deriv_index,"*",err,error))// &
3845  & " is invalid."
3846  CALL flagerror(local_error,err,error,*999)
3847  END SELECT
3848  CASE(field_deludeln_variable_type)
3849  SELECT CASE(global_deriv_index)
3850  CASE(no_global_deriv)
3851  IF(component_idx==1) THEN
3852  !calculate u
3853  VALUE=0.0_dp
3854  ELSE IF(component_idx==2) THEN
3855  !calculate v
3856  VALUE=16.0_dp*mu_param*pi**2/10.0_dp**2*cos(2.0_dp*pi*x(2)/10.0_dp)*cos(2.0_dp*pi*x(1)/10.0_dp)
3857  ELSE IF(component_idx==3) THEN
3858  !calculate p
3859  VALUE=0.0_dp
3860  ELSE
3861  CALL flagerror("Not implemented.",err,error,*999)
3862  ENDIF
3863  CASE(global_deriv_s1)
3864  CALL flagerror("Not implemented.",err,error,*999)
3865  CASE(global_deriv_s2)
3866  CALL flagerror("Not implemented.",err,error,*999)
3867  CASE(global_deriv_s1_s2)
3868  CALL flagerror("Not implemented.",err,error,*999)
3869  CASE DEFAULT
3870  local_error="The global derivative index of "//trim(number_to_vstring( &
3871  & global_deriv_index,"*",err,error))// &
3872  & " is invalid."
3873  CALL flagerror(local_error,err,error,*999)
3874  END SELECT
3875  CASE DEFAULT
3876  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
3877  & " is invalid."
3878  CALL flagerror(local_error,err,error,*999)
3879  END SELECT
3880  ELSE
3881  local_error="The number of components does not correspond to the number of dimensions."
3882  CALL flagerror(local_error,err,error,*999)
3883  ENDIF
3884  CASE(equations_set_stokes_equation_two_dim_4)
3885  IF(number_of_dimensions==2.AND.number_of_components==3) THEN
3886  !Reduced Taylor-Green solution for Stokes
3887  CALL flagerror("Not implemented.",err,error,*999)
3888  ENDIF
3889  CASE(equations_set_stokes_equation_two_dim_5)
3890  IF(number_of_dimensions==2.AND.number_of_components==3) THEN
3891  !Stokes-Taylor-Green dynamic
3892  SELECT CASE(variable_type)
3893  CASE(field_u_variable_type)
3894  SELECT CASE(global_deriv_index)
3895  CASE(no_global_deriv)
3896  IF(component_idx==1) THEN
3897  !calculate u
3898  VALUE=x(2)*exp(-(2.0_dp*mu_param/rho_param*current_time))
3899  ELSE IF(component_idx==2) THEN
3900  !calculate v
3901  VALUE=x(1)*exp(-(2.0_dp*mu_param/rho_param*current_time))
3902  ELSE IF(component_idx==3) THEN
3903  !calculate p
3904  VALUE=2.0_dp*x(2)*mu_param*exp(-(2.0_dp*mu_param/rho_param*current_time))*x(1)
3905  ELSE
3906  CALL flagerror("Not implemented.",err,error,*999)
3907  ENDIF
3908  CASE(global_deriv_s1)
3909  CALL flagerror("Not implemented.",err,error,*999)
3910  CASE(global_deriv_s2)
3911  CALL flagerror("Not implemented.",err,error,*999)
3912  CASE(global_deriv_s1_s2)
3913  CALL flagerror("Not implemented.",err,error,*999)
3914  CASE DEFAULT
3915  local_error="The global derivative index of "//trim(number_to_vstring( &
3916  & global_deriv_index,"*",err,error))// &
3917  & " is invalid."
3918  CALL flagerror(local_error,err,error,*999)
3919  END SELECT
3920  CASE(field_deludeln_variable_type)
3921  SELECT CASE(global_deriv_index)
3922  CASE(no_global_deriv)
3923  IF(component_idx==1) THEN
3924  !calculate u
3925  VALUE=0.0_dp
3926  ELSE IF(component_idx==2) THEN
3927  !calculate v
3928  VALUE=0.0_dp
3929  ELSE IF(component_idx==3) THEN
3930  !calculate p
3931  VALUE=0.0_dp
3932  ELSE
3933  CALL flagerror("Not implemented.",err,error,*999)
3934  ENDIF
3935  CASE(global_deriv_s1)
3936  CALL flagerror("Not implemented.",err,error,*999)
3937  CASE(global_deriv_s2)
3938  CALL flagerror("Not implemented.",err,error,*999)
3939  CASE(global_deriv_s1_s2)
3940  CALL flagerror("Not implemented.",err,error,*999)
3941  CASE DEFAULT
3942  local_error="The global derivative index of "//trim(number_to_vstring( &
3943  & global_deriv_index,"*",err,error))// &
3944  & " is invalid."
3945  CALL flagerror(local_error,err,error,*999)
3946  END SELECT
3947  CASE DEFAULT
3948  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
3949  & " is invalid."
3950  CALL flagerror(local_error,err,error,*999)
3951  END SELECT
3952  ELSE
3953  local_error="The number of components does not correspond to the number of dimensions."
3954  CALL flagerror(local_error,err,error,*999)
3955  ENDIF
3956  CASE(equations_set_stokes_equation_three_dim_1)
3957  IF(number_of_dimensions==3.AND.number_of_components==4) THEN
3958 !POLYNOM
3959  SELECT CASE(variable_type)
3960  CASE(field_u_variable_type)
3961  SELECT CASE(global_deriv_index)
3962  CASE(no_global_deriv)
3963  IF(component_idx==1) THEN
3964  !calculate u
3965  VALUE=x(2)**2/10.0_dp**2+x(3)**2/10.0_dp**2
3966  ELSE IF(component_idx==2) THEN
3967  !calculate v
3968  VALUE=x(1)**2/10.0_dp**2+x(3)**2/10.0_dp**2
3969  ELSE IF(component_idx==3) THEN
3970  !calculate w
3971  VALUE=x(1)**2/10.0_dp**2+x(2)**2/10.0_dp**2
3972  ELSE IF(component_idx==4) THEN
3973  !calculate p
3974  VALUE=4.0_dp*mu_param/10.0_dp**2*x(1)
3975  ELSE
3976  CALL flagerror("Not implemented.",err,error,*999)
3977  ENDIF
3978  CASE(global_deriv_s1)
3979  CALL flagerror("Not implemented.",err,error,*999)
3980  CASE(global_deriv_s2)
3981  CALL flagerror("Not implemented.",err,error,*999)
3982  CASE(global_deriv_s1_s2)
3983  CALL flagerror("Not implemented.",err,error,*999)
3984  CASE DEFAULT
3985  local_error="The global derivative index of "//trim(number_to_vstring( &
3986  & global_deriv_index,"*",err,error))// &
3987  & " is invalid."
3988  CALL flagerror(local_error,err,error,*999)
3989  END SELECT
3990  CASE(field_deludeln_variable_type)
3991  SELECT CASE(global_deriv_index)
3992  CASE(no_global_deriv)
3993  VALUE=0.0_dp
3994  CASE(global_deriv_s1)
3995  CALL flagerror("Not implemented.",err,error,*999)
3996  CASE(global_deriv_s2)
3997  CALL flagerror("Not implemented.",err,error,*999)
3998  CASE(global_deriv_s1_s2)
3999  CALL flagerror("Not implemented.",err,error,*999)
4000  CASE DEFAULT
4001  local_error="The global derivative index of "//trim(number_to_vstring( &
4002  & global_deriv_index,"*",err,error))// &
4003  & " is invalid."
4004  CALL flagerror(local_error,err,error,*999)
4005  END SELECT
4006  CASE DEFAULT
4007  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
4008  & " is invalid."
4009  CALL flagerror(local_error,err,error,*999)
4010  END SELECT
4011  ELSE
4012  local_error="The number of components does not correspond to the number of dimensions."
4013  CALL flagerror(local_error,err,error,*999)
4014  ENDIF
4015  CASE(equations_set_stokes_equation_three_dim_2)
4016  IF(number_of_dimensions==3.AND.number_of_components==4) THEN
4017  !Exponential function
4018  SELECT CASE(variable_type)
4019  CASE(field_u_variable_type)
4020  SELECT CASE(global_deriv_index)
4021  CASE(no_global_deriv)
4022  IF(component_idx==1) THEN
4023  !calculate u
4024  VALUE=exp((x(1)-x(2))/10.0_dp)+exp((x(3)-x(1))/10.0_dp)
4025  ELSE IF(component_idx==2) THEN
4026  !calculate v
4027  VALUE=exp((x(1)-x(2))/10.0_dp)+exp((x(2)-x(3))/10.0_dp)
4028  ELSE IF(component_idx==3) THEN
4029  !calculate w
4030  VALUE=exp((x(3)-x(1))/10.0_dp)+exp((x(2)-x(3))/10.0_dp)
4031  ELSE IF(component_idx==4) THEN
4032  !calculate p
4033  VALUE=2.0_dp*mu_param/10.0_dp*(exp((x(1)-x(2))/10.0_dp)-exp((x(3)-x(1))/10.0_dp))
4034  ELSE
4035  CALL flagerror("Not implemented.",err,error,*999)
4036  ENDIF
4037  CASE(global_deriv_s1)
4038  CALL flagerror("Not implemented.",err,error,*999)
4039  CASE(global_deriv_s2)
4040  CALL flagerror("Not implemented.",err,error,*999)
4041  CASE(global_deriv_s1_s2)
4042  CALL flagerror("Not implemented.",err,error,*999)
4043  CASE DEFAULT
4044  local_error="The global derivative index of "//trim(number_to_vstring( &
4045  & global_deriv_index,"*",err,error))// &
4046  & " is invalid."
4047  CALL flagerror(local_error,err,error,*999)
4048  END SELECT
4049  CASE(field_deludeln_variable_type)
4050  SELECT CASE(global_deriv_index)
4051  CASE(no_global_deriv)
4052  IF(component_idx==1) THEN
4053  !calculate u
4054  VALUE=0.0_dp
4055  ELSE IF(component_idx==2) THEN
4056  !calculate v
4057  VALUE=-2.0_dp*mu_param*(2.0_dp*exp(x(1)-x(2))+exp(x(2)-x(3)))
4058  ELSE IF(component_idx==3) THEN
4059  !calculate w
4060  VALUE=-2.0_dp*mu_param*(2.0_dp*exp(x(3)-x(1))+exp(x(2)-x(3)))
4061  ELSE IF(component_idx==4) THEN
4062  !calculate p
4063  VALUE=0.0_dp
4064  ELSE
4065  CALL flagerror("Not implemented.",err,error,*999)
4066  ENDIF
4067  CASE(global_deriv_s1)
4068  CALL flagerror("Not implemented.",err,error,*999)
4069  CASE(global_deriv_s2)
4070  CALL flagerror("Not implemented.",err,error,*999)
4071  CASE(global_deriv_s1_s2)
4072  CALL flagerror("Not implemented.",err,error,*999)
4073  CASE DEFAULT
4074  local_error="The global derivative index of "//trim(number_to_vstring( &
4075  & global_deriv_index,"*",err,error))// &
4076  & " is invalid."
4077  CALL flagerror(local_error,err,error,*999)
4078  END SELECT
4079  CASE DEFAULT
4080  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
4081  & " is invalid."
4082  CALL flagerror(local_error,err,error,*999)
4083  END SELECT
4084  ELSE
4085  local_error="The number of components does not correspond to the number of dimensions."
4086  CALL flagerror(local_error,err,error,*999)
4087  ENDIF
4088  CASE(equations_set_stokes_equation_three_dim_3)
4089  IF(number_of_dimensions==3.AND.number_of_components==4) THEN
4090  !Sine and cosine functions
4091  SELECT CASE(variable_type)
4092  CASE(field_u_variable_type)
4093  SELECT CASE(global_deriv_index)
4094  CASE(no_global_deriv)
4095  IF(component_idx==1) THEN
4096  !calculate u
4097  VALUE=sin(2.0_dp*pi*x(1)/10.0_dp)*sin(2.0_dp*pi*x(2)/10.0_dp)*sin(2.0_dp*pi*x(3)/10.0_dp)
4098  ELSE IF(component_idx==2) THEN
4099  !calculate v
4100  VALUE=2.0_dp*cos(2.0_dp*pi*x(1)/10.0_dp)*sin(2.0_dp*pi*x(3)/10.0_dp)*cos(2.0_dp*pi*x(2)/10.0_dp)
4101  ELSE IF(component_idx==3) THEN
4102  !calculate w
4103  VALUE=-cos(2.0_dp*pi*x(1)/10.0_dp)*sin(2.0_dp*pi*x(2)/10.0_dp)*cos(2.0_dp*pi*x(3)/10.0_dp)
4104  ELSE IF(component_idx==4) THEN
4105  !calculate p
4106  VALUE=6.0_dp*mu_param*pi/10.0_dp*sin(2.0_dp*pi*x(2)/10.0_dp)*sin(2.0_dp*pi*x(3)/10.0_dp)* &
4107  & cos(2.0_dp*pi*x(1)/10.0_dp)
4108  ELSE
4109  CALL flagerror("Not implemented.",err,error,*999)
4110  ENDIF
4111  CASE(global_deriv_s1)
4112  CALL flagerror("Not implemented.",err,error,*999)
4113  CASE(global_deriv_s2)
4114  CALL flagerror("Not implemented.",err,error,*999)
4115  CASE(global_deriv_s1_s2)
4116  CALL flagerror("Not implemented.",err,error,*999)
4117  CASE DEFAULT
4118  local_error="The global derivative index of "//trim(number_to_vstring( &
4119  & global_deriv_index,"*",err,error))// &
4120  & " is invalid."
4121  CALL flagerror(local_error,err,error,*999)
4122  END SELECT
4123  CASE(field_deludeln_variable_type)
4124  SELECT CASE(global_deriv_index)
4125  CASE(no_global_deriv)
4126  IF(component_idx==1) THEN
4127  !calculate u
4128  VALUE=0.0_dp
4129  ELSE IF(component_idx==2) THEN
4130  !calculate v
4131  VALUE=36*mu_param*pi**2/10.0_dp**2*cos(2.0_dp*pi*x(2)/10.0_dp)*sin(2.0_dp*pi*x(3)/10.0_dp)* &
4132  & cos(2.0_dp*pi*x(1)/10.0_dp)
4133  ELSE IF(component_idx==3) THEN
4134  !calculate w
4135  VALUE=0.0_dp
4136  ELSE IF(component_idx==4) THEN
4137  !calculate p
4138  VALUE=0.0_dp
4139  ELSE
4140  CALL flagerror("Not implemented.",err,error,*999)
4141  ENDIF
4142  CASE(global_deriv_s1)
4143  CALL flagerror("Not implemented.",err,error,*999)
4144  CASE(global_deriv_s2)
4145  CALL flagerror("Not implemented.",err,error,*999)
4146  CASE(global_deriv_s1_s2)
4147  CALL flagerror("Not implemented.",err,error,*999)
4148  CASE DEFAULT
4149  local_error="The global derivative index of "//trim(number_to_vstring( &
4150  & global_deriv_index,"*",err,error))// &
4151  & " is invalid."
4152  CALL flagerror(local_error,err,error,*999)
4153  END SELECT
4154  CASE DEFAULT
4155  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
4156  & " is invalid."
4157  CALL flagerror(local_error,err,error,*999)
4158  END SELECT
4159  ELSE
4160  local_error="The number of components does not correspond to the number of dimensions."
4161  CALL flagerror(local_error,err,error,*999)
4162  ENDIF
4163  CASE(equations_set_stokes_equation_three_dim_4)
4164  IF(number_of_dimensions==3.AND.number_of_components==4) THEN
4165  !Reduced Taylor-Green solution for Stokes
4166  CALL flagerror("Not implemented.",err,error,*999)
4167  ENDIF
4168  CASE(equations_set_stokes_equation_three_dim_5)
4169  IF(number_of_dimensions==3.AND.number_of_components==4) THEN
4170  !Stokes-Taylor-Green dynamic
4171  SELECT CASE(variable_type)
4172  CASE(field_u_variable_type)
4173  SELECT CASE(global_deriv_index)
4174  CASE(no_global_deriv)
4175  IF(component_idx==1) THEN
4176  !calculate u
4177  VALUE=x(2)*exp(-(2.0_dp*mu_param/rho_param*current_time))
4178  ELSE IF(component_idx==2) THEN
4179  !calculate v
4180  VALUE=x(1)*exp(-(2.0_dp*mu_param/rho_param*current_time))
4181  ELSE IF(component_idx==3) THEN
4182  !calculate v
4183  VALUE=0.0_dp
4184  ELSE IF(component_idx==4) THEN
4185  !calculate p
4186  VALUE=2.0_dp*x(2)*mu_param*exp(-(2.0_dp*mu_param/rho_param*current_time))*x(1)
4187  ELSE
4188  CALL flagerror("Not implemented.",err,error,*999)
4189  ENDIF
4190  CASE(global_deriv_s1)
4191  CALL flagerror("Not implemented.",err,error,*999)
4192  CASE(global_deriv_s2)
4193  CALL flagerror("Not implemented.",err,error,*999)
4194  CASE(global_deriv_s1_s2)
4195  CALL flagerror("Not implemented.",err,error,*999)
4196  CASE DEFAULT
4197  local_error="The global derivative index of "//trim(number_to_vstring( &
4198  & global_deriv_index,"*",err,error))// &
4199  & " is invalid."
4200  CALL flagerror(local_error,err,error,*999)
4201  END SELECT
4202  CASE(field_deludeln_variable_type)
4203  SELECT CASE(global_deriv_index)
4204  CASE(no_global_deriv)
4205  IF(component_idx==1) THEN
4206  !calculate u
4207  VALUE=0.0_dp
4208  ELSE IF(component_idx==2) THEN
4209  !calculate v
4210  VALUE=0.0_dp
4211  ELSE IF(component_idx==3) THEN
4212  !calculate p
4213  VALUE=0.0_dp
4214  ELSE IF(component_idx==4) THEN
4215  !calculate p
4216  VALUE=0.0_dp
4217  ELSE
4218  CALL flagerror("Not implemented.",err,error,*999)
4219  ENDIF
4220  CASE(global_deriv_s1)
4221  CALL flagerror("Not implemented.",err,error,*999)
4222  CASE(global_deriv_s2)
4223  CALL flagerror("Not implemented.",err,error,*999)
4224  CASE(global_deriv_s1_s2)
4225  CALL flagerror("Not implemented.",err,error,*999)
4226  CASE DEFAULT
4227  local_error="The global derivative index of "//trim(number_to_vstring( &
4228  & global_deriv_index,"*",err,error))// &
4229  & " is invalid."
4230  CALL flagerror(local_error,err,error,*999)
4231  END SELECT
4232  CASE DEFAULT
4233  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
4234  & " is invalid."
4235  CALL flagerror(local_error,err,error,*999)
4236  END SELECT
4237  ELSE
4238  local_error="The number of components does not correspond to the number of dimensions."
4239  CALL flagerror(local_error,err,error,*999)
4240  ENDIF
4241  CASE DEFAULT
4242  local_error="The analytic function type of "// &
4243  & trim(number_to_vstring(analytic_function_type,"*",err,error))// &
4244  & " is invalid."
4245  CALL flagerror(local_error,err,error,*999)
4246  END SELECT
4247  exits("STOKES_EQUATION_ANALYTIC_FUNCTIONS")
4248  RETURN
4249 999 errorsexits("STOKES_EQUATION_ANALYTIC_FUNCTIONS",err,error)
4250  RETURN 1
4251  END SUBROUTINE stokes_equation_analytic_functions
4252  !
4253  !================================================================================================================================
4254  !
4255 
4256 END MODULE stokes_equations_routines
4257 
4258  !
4259  !================================================================================================================================
4260  !
equations_set_constants::equations_set_setup_dependent_type
integer(intg), parameter equations_set_setup_dependent_type
Dependent variables.
Definition: equations_set_constants.f90:324
equations_set_constants::equations_set_fem_solution_method
integer(intg), parameter equations_set_fem_solution_method
Finite Element Method solution method.
Definition: equations_set_constants.f90:348
basis_routines
This module contains all basis function routines.
Definition: basis_routines.f90:45
equations_set_constants::equations_set_setup_materials_type
integer(intg), parameter equations_set_setup_materials_type
Materials setup.
Definition: equations_set_constants.f90:326
types::boundary_conditions_type
Contains information on the boundary conditions for the solver equations.
Definition: types.f90:1780
base_routines::enters
subroutine, public enters(NAME, ERR, ERROR,)
Records the entry into the named procedure and initialises the error code.
Definition: base_routines.f90:363
boundary_conditions_routines::boundary_condition_moved_wall
integer(intg), parameter, public boundary_condition_moved_wall
The dof is fixed as a boundary condition.
Definition: boundary_condition_routines.f90:99
solver_routines::solvers_create_finish
subroutine, public solvers_create_finish(SOLVERS, ERR, ERROR,)
Finish the creation of solvers.
Definition: solver_routines.f90:22598
kinds::ptr
integer, parameter ptr
Pointer integer kind.
Definition: kinds.f90:58
control_loop_routines::control_loop_progress_output
integer(intg), parameter, public control_loop_progress_output
Progress output from control loop.
Definition: control_loop_routines.f90:75
equations_mapping_routines::equations_mapping_dynamic_variable_type_set
subroutine, public equations_mapping_dynamic_variable_type_set(EQUATIONS_MAPPING, DYNAMIC_VARIABLE_TYPE, ERR, ERROR,)
Sets the mapping between a dependent field variable and the equations set dynamic matrices...
Definition: equations_mapping_routines.f90:1996
types::equations_mapping_type
Contains information on the equations mapping i.e., how field variable DOFS are mapped to the rows an...
Definition: types.f90:1681
types::equations_type
Contains information about the equations in an equations set.
Definition: types.f90:1735
equations_set_constants::equations_set_gfem_solution_method
integer(intg), parameter equations_set_gfem_solution_method
Grid-based Finite Element Method solution method.
Definition: equations_set_constants.f90:352
problem_constants::problem_control_time_loop_type
integer(intg), parameter problem_control_time_loop_type
Time control loop.
Definition: problem_constants.f90:283
problem_constants::problem_setup_control_type
integer(intg), parameter problem_setup_control_type
Solver setup for a problem.
Definition: problem_constants.f90:263
problem_constants
This module handles all problem wide constants.
Definition: problem_constants.f90:45
problem_constants::solver_equations_first_order_dynamic
integer(intg), parameter solver_equations_first_order_dynamic
Solver equations are first order dynamic.
Definition: problem_constants.f90:311
control_loop_routines::control_loop_node
integer(intg), parameter, public control_loop_node
The identifier for a each "leaf" node in a control loop.
Definition: control_loop_routines.f90:67
solver_routines::solver_dynamic_order_set
subroutine, public solver_dynamic_order_set(SOLVER, ORDER, ERR, ERROR,)
Sets/changes the order for a dynamic solver.
Definition: solver_routines.f90:6327
constants::no_global_deriv
integer(intg), parameter no_global_deriv
No global derivative i.e., u.
Definition: constants.f90:213
strings::number_to_vstring
Converts a number to its equivalent varying string representation.
Definition: strings.f90:161
equations_routines::equations_create_start
subroutine, public equations_create_start(EQUATIONS_SET, EQUATIONS, ERR, ERROR,)
Start the creation of equations for the equation set.
Definition: equations_routines.f90:192
types::decomposition_type
Contains information on the mesh decomposition.
Definition: types.f90:1063
equations_set_constants::equations_set_transient_stokes_subtype
integer(intg), parameter equations_set_transient_stokes_subtype
Definition: equations_set_constants.f90:162
equations_set_constants::equations_set_stokes_equation_two_dim_3
integer(intg), parameter equations_set_stokes_equation_two_dim_3
u=tbd
Definition: equations_set_constants.f90:490
equations_matrices_routines::equations_matrices_create_start
subroutine, public equations_matrices_create_start(EQUATIONS, EQUATIONS_MATRICES, ERR, ERROR,)
Starts the creation of the equations matrices and rhs for the the equations.
Definition: equations_matrices_routines.f90:615
types::solver_type
Contains information on the type of solver to be used.
Definition: types.f90:2777
solver_routines::solver_petsc_library
integer(intg), parameter, public solver_petsc_library
PETSc solver library.
Definition: solver_routines.f90:104
constants::pi
real(dp), parameter pi
The double precision value of pi.
Definition: constants.f90:57
solver_routines::solvers_number_set
subroutine, public solvers_number_set(SOLVERS, NUMBER_OF_SOLVERS, ERR, ERROR,)
Sets/changes the number of solvers.
Definition: solver_routines.f90:22805
constants::no_part_deriv
integer(intg), parameter no_part_deriv
No partial derivative i.e., u.
Definition: constants.f90:177
solver_routines::solver_dynamic_crank_nicolson_scheme
integer(intg), parameter, public solver_dynamic_crank_nicolson_scheme
Crank-Nicolson dynamic solver.
Definition: solver_routines.f90:282
solver_routines::solver_dynamic_degree_set
subroutine, public solver_dynamic_degree_set(SOLVER, DEGREE, ERR, ERROR,)
Sets/changes the degree of the polynomial used to interpolate time for a dynamic solver.
Definition: solver_routines.f90:5484
stokes_equations_routines::stokes_problem_setup
subroutine, public stokes_problem_setup(PROBLEM, PROBLEM_SETUP, ERR, ERROR,)
Sets up the Stokes problem.
Definition: Stokes_equations_routines.f90:981
stokes_equations_routines::stokes_equations_set_setup
subroutine, public stokes_equations_set_setup(EQUATIONS_SET, EQUATIONS_SET_SETUP, ERR, ERROR,)
Sets up the standard Stokes fluid setup.
Definition: Stokes_equations_routines.f90:230
equations_matrices_routines
This module handles all equations matrix and rhs routines.
Definition: equations_matrices_routines.f90:45
solver_routines::solver_dynamic_first_order
integer(intg), parameter, public solver_dynamic_first_order
Dynamic solver has first order terms.
Definition: solver_routines.f90:255
equations_set_constants::equations_set_navier_stokes_equation_two_dim_4
integer(intg), parameter equations_set_navier_stokes_equation_two_dim_4
u=tbd
Definition: equations_set_constants.f90:505
stokes_equations_routines::stokes_pre_solve
subroutine, public stokes_pre_solve(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Sets up the Stokes problem pre solve.
Definition: Stokes_equations_routines.f90:2004
solver_routines::solver_type_set
subroutine, public solver_type_set(SOLVER, SOLVE_TYPE, ERR, ERROR,)
Sets/changes the type for a solver.
Definition: solver_routines.f90:21345
stokes_equations_routines::stokes_equation_analytic_functions
subroutine, public stokes_equation_analytic_functions(VALUE, X, MU_PARAM, RHO_PARAM, CURRENT_TIME, VARIABLE_TYPE, GLOBAL_DERIV_INDEX, ANALYTIC_FUNCTION_TYPE, NUMBER_OF_DIMENSIONS, NUMBER_OF_COMPONENTS, COMPONENT_IDX, ERR, ERROR,)
Calculates the various analytic solutions given X and time, can be called from within analytic calcul...
Definition: Stokes_equations_routines.f90:3673
equations_set_constants::equations_static
integer(intg), parameter equations_static
The equations are static and have no time dependence.
Definition: equations_set_constants.f90:376
types::equations_set_type
Contains information on an equations set.
Definition: types.f90:1941
equations_routines
This module handles all equations routines.
Definition: equations_routines.f90:45
equations_set_constants::equations_set_setup_source_type
integer(intg), parameter equations_set_setup_source_type
Source setup.
Definition: equations_set_constants.f90:327
types::fields_type
Contains information on the fields defined on a region.
Definition: types.f90:1373
strings
This module contains all string manipulation and transformation routines.
Definition: strings.f90:45
solver_routines::solvers_create_start
subroutine, public solvers_create_start(CONTROL_LOOP, SOLVERS, ERR, ERROR,)
Start the creation of a solvers for the control loop.
Definition: solver_routines.f90:22651
types::solvers_type
Contains information on the solvers to be used in a control loop.
Definition: types.f90:2805
equations_mapping_routines::equations_mapping_dynamic_matrices_set
Definition: equations_mapping_routines.f90:83
constants::first_part_deriv
integer(intg), parameter first_part_deriv
First partial derivative i.e., du/ds.
Definition: constants.f90:178
timer
This module contains routines for timing the program.
Definition: timer_f.f90:45
control_loop_routines::control_loop_current_times_get
subroutine, public control_loop_current_times_get(CONTROL_LOOP, CURRENT_TIME, TIME_INCREMENT, ERR, ERROR,)
Gets the current time parameters for a time control loop.
Definition: control_loop_routines.f90:231
problem_constants::problem_ale_stokes_subtype
integer(intg), parameter problem_ale_stokes_subtype
Definition: problem_constants.f90:136
equations_matrices_routines::equations_matrices_dynamic_lumping_type_set
subroutine, public equations_matrices_dynamic_lumping_type_set(EQUATIONS_MATRICES, LUMPING_TYPE, ERR, ERROR,)
Sets the lumping of the linear equations matrices.
Definition: equations_matrices_routines.f90:3760
problem_constants::solver_equations_static
integer(intg), parameter solver_equations_static
Solver equations are static.
Definition: problem_constants.f90:309
equations_routines::equations_time_dependence_type_set
subroutine, public equations_time_dependence_type_set(EQUATIONS, TIME_DEPENDENCE_TYPE, ERR, ERROR,)
Sets/changes the time dependence type for equations.
Definition: equations_routines.f90:851
problem_constants::problem_stokes_equation_type
integer(intg), parameter problem_stokes_equation_type
Definition: problem_constants.f90:73
solver_routines::solver_equations_sparsity_type_set
subroutine, public solver_equations_sparsity_type_set(SOLVER_EQUATIONS, SPARSITY_TYPE, ERR, ERROR,)
Sets/changes the sparsity type for solver equations.
Definition: solver_routines.f90:7734
analytic_analysis_routines
This module handles all analytic analysis routines.
Definition: analytic_analysis_routines.f90:45
solver_routines::solvers_solver_get
subroutine, public solvers_solver_get(SOLVERS, SOLVER_INDEX, SOLVER, ERR, ERROR,)
Returns a pointer to the specified solver in the list of solvers.
Definition: solver_routines.f90:22876
stokes_equations_routines::stokes_equationssetsolutionmethodset
subroutine, public stokes_equationssetsolutionmethodset(EQUATIONS_SET, SOLUTION_METHOD, ERR, ERROR,)
Sets/changes the solution method for a Stokes flow equation type of an fluid mechanics equations set ...
Definition: Stokes_equations_routines.f90:108
types::field_type
Contains information for a field defined on a region.
Definition: types.f90:1346
stokes_equations_routines::stokes_pre_solve_ale_update_mesh
subroutine stokes_pre_solve_ale_update_mesh(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Update mesh velocity and move mesh for ALE Stokes problem.
Definition: Stokes_equations_routines.f90:2677
equations_matrices_routines::equations_matrices_full_matrices
integer(intg), parameter, public equations_matrices_full_matrices
Use fully populated equation matrices.
Definition: equations_matrices_routines.f90:91
equations_set_constants::equations_set_navier_stokes_equation_two_dim_5
integer(intg), parameter equations_set_navier_stokes_equation_two_dim_5
u=tbd
Definition: equations_set_constants.f90:506
equations_set_constants::equations_set_fluid_mechanics_class
integer(intg), parameter equations_set_fluid_mechanics_class
Definition: equations_set_constants.f90:55
equations_mapping_routines::equations_mapping_rhs_variable_type_set
subroutine, public equations_mapping_rhs_variable_type_set(EQUATIONS_MAPPING, RHS_VARIABLE_TYPE, ERR, ERROR,)
Sets the mapping between a dependent field variable and the equations set rhs vector.
Definition: equations_mapping_routines.f90:3170
problem_constants::solver_equations_linear
integer(intg), parameter solver_equations_linear
Solver equations are linear.
Definition: problem_constants.f90:301
constants::global_deriv_s2
integer(intg), parameter global_deriv_s2
First global derivative in the s2 direction i.e., du/ds2.
Definition: constants.f90:215
types::control_loop_type
Contains information on a control loop.
Definition: types.f90:3185
iso_varying_string
This module provides an iso_varying_string module, conformant to the API specified in ISO/IEC 1539-2:...
Definition: iso_varying_string.f90:62
solver_routines::solver_equations_create_finish
subroutine, public solver_equations_create_finish(SOLVER_EQUATIONS, ERR, ERROR,)
Finishes the process of creating solver equations.
Definition: solver_routines.f90:7166
solver_routines::solver_sparse_matrices
integer(intg), parameter, public solver_sparse_matrices
Use sparse solver matrices.
Definition: solver_routines.f90:360
stokes_equations_routines::stokes_post_solve
subroutine, public stokes_post_solve(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Sets up the Stokes problem post solve.
Definition: Stokes_equations_routines.f90:1934
solver_routines::solver_equations_create_start
subroutine, public solver_equations_create_start(SOLVER, SOLVER_EQUATIONS, ERR, ERROR,)
Starts the process of creating solver equations.
Definition: solver_routines.f90:7209
stokes_equations_routines::stokes_pre_solve_update_boundary_conditions
subroutine stokes_pre_solve_update_boundary_conditions(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Update boundary conditions for Stokes flow pre solve.
Definition: Stokes_equations_routines.f90:2092
solver_routines::solver_dynamic_type
integer(intg), parameter, public solver_dynamic_type
A dynamic solver.
Definition: solver_routines.f90:90
equations_set_constants::equations_set_stokes_equation_three_dim_4
integer(intg), parameter equations_set_stokes_equation_three_dim_4
u=tbd
Definition: equations_set_constants.f90:496
basis_routines::basis_default_quadrature_scheme
integer(intg), parameter, public basis_default_quadrature_scheme
Identifier for the default quadrature scheme.
Definition: basis_routines.f90:125
problem_constants::problem_setup_solvers_type
integer(intg), parameter problem_setup_solvers_type
Solver setup for a problem.
Definition: problem_constants.f90:264
equations_set_constants::equations_set_setup_equations_type
integer(intg), parameter equations_set_setup_equations_type
Equations setup.
Definition: equations_set_constants.f90:330
types::equations_mapping_dynamic_type
Contains information for mapping field variables to the dynamic matrices in the equations set of the ...
Definition: types.f90:1571
iso_varying_string::varying_string
Definition: iso_varying_string.f90:75
equations_set_constants::equations_set_setup_independent_type
integer(intg), parameter equations_set_setup_independent_type
Independent variables.
Definition: equations_set_constants.f90:325
constants
This module contains all program wide constants.
Definition: constants.f90:45
solver_routines::solver_library_type_set
subroutine, public solver_library_type_set(SOLVER, SOLVER_LIBRARY_TYPE, ERR, ERROR,)
Sets/changes the type of library type to use for the solver.
Definition: solver_routines.f90:9067
boundary_conditions_routines::boundary_condition_fixed_inlet
integer(intg), parameter, public boundary_condition_fixed_inlet
The dof is fixed as a boundary condition.
Definition: boundary_condition_routines.f90:96
equations_mapping_routines::equationsmapping_linearmatricesnumberset
subroutine, public equationsmapping_linearmatricesnumberset(EQUATIONS_MAPPING, NUMBER_OF_LINEAR_EQUATIONS_MATRICES, ERR, ERROR,)
Sets/changes the number of linear equations matrices.
Definition: equations_mapping_routines.f90:2484
equations_routines::equations_lumped_matrices
integer(intg), parameter, public equations_lumped_matrices
The equations matrices are "mass" lumped.
Definition: equations_routines.f90:88
equations_set_constants::equations_set_laplace_stokes_subtype
integer(intg), parameter equations_set_laplace_stokes_subtype
Definition: equations_set_constants.f90:161
equations_set_constants::equations_set_stokes_equation_two_dim_2
integer(intg), parameter equations_set_stokes_equation_two_dim_2
u=tbd
Definition: equations_set_constants.f90:489
equations_set_constants::equations_set_navier_stokes_equation_three_dim_5
integer(intg), parameter equations_set_navier_stokes_equation_three_dim_5
u=tbd
Definition: equations_set_constants.f90:511
problem_constants::problem_setup_initial_type
integer(intg), parameter problem_setup_initial_type
Initial setup for a problem.
Definition: problem_constants.f90:262
equations_mapping_routines::equationsmapping_linearmatricesvariabletypesset
subroutine, public equationsmapping_linearmatricesvariabletypesset(EQUATIONS_MAPPING, LINEAR_MATRIX_VARIABLE_TYPES, ERR, ERROR,)
Sets the mapping between the dependent field variable types and the linear equations matrices...
Definition: equations_mapping_routines.f90:2716
types::problem_setup_type
Definition: types.f90:3215
fluid_mechanics_io_routines::fluid_mechanics_io_read_data
subroutine, public fluid_mechanics_io_read_data(SOLVER_TYPE, INPUT_VALUES, NUMBER_OF_DIMENSIONS, INPUT_TYPE, INPUT_OPTION, TIME_STEP, LENGTH_SCALE)
Reads input data from a file.
Definition: fluid_mechanics_IO_routines.f90:4319
equations_set_constants::equations_set_stokes_equation_two_dim_4
integer(intg), parameter equations_set_stokes_equation_two_dim_4
u=tbd
Definition: equations_set_constants.f90:491
equations_set_constants::equations_first_order_dynamic
integer(intg), parameter equations_first_order_dynamic
The equations are first order dynamic.
Definition: equations_set_constants.f90:378
types::boundary_conditions_variable_type
Contains information on the boundary conditions for a dependent field variable.
Definition: types.f90:1759
equations_set_constants::equations_set_static_stokes_subtype
integer(intg), parameter equations_set_static_stokes_subtype
Definition: equations_set_constants.f90:160
solver_routines::solver_equations_linearity_type_set
subroutine, public solver_equations_linearity_type_set(SOLVER_EQUATIONS, LINEARITY_TYPE, ERR, ERROR,)
Sets/changes the linearity type for solver equations.
Definition: solver_routines.f90:7577
problem_constants::problem_laplace_stokes_subtype
integer(intg), parameter problem_laplace_stokes_subtype
Definition: problem_constants.f90:133
equations_set_constants::equations_set_setup_start_action
integer(intg), parameter equations_set_setup_start_action
Start setup action.
Definition: equations_set_constants.f90:338
base_routines::exits
subroutine, public exits(NAME)
Records the exit out of the named procedure.
Definition: base_routines.f90:504
control_loop_routines::control_loop_solvers_get
recursive subroutine, public control_loop_solvers_get(CONTROL_LOOP, SOLVERS, ERR, ERROR,)
Returns a pointer to the solvers for a control loop.
Definition: control_loop_routines.f90:1333
types
This module contains all type definitions in order to avoid cyclic module references.
Definition: types.f90:70
types::equations_matrices_type
Contains information on the equations matrices and vectors.
Definition: types.f90:1520
equations_matrices_routines::equations_matrix_fem_structure
integer(intg), parameter, public equations_matrix_fem_structure
Finite element matrix structure.
Definition: equations_matrices_routines.f90:72
input_output::write_string
Write a string to a given output stream.
Definition: input_output.f90:74
base_routines
This module contains all the low-level base routines e.g., all debug, control, and low-level communic...
Definition: base_routines.f90:45
types::equations_matrices_linear_type
Contains information of the linear matrices for equations matrices.
Definition: types.f90:1479
base_routines::general_output_type
integer(intg), parameter, public general_output_type
General output type.
Definition: base_routines.f90:66
equations_matrices_routines::equationsmatrices_dynamicstructuretypeset
subroutine, public equationsmatrices_dynamicstructuretypeset(EQUATIONS_MATRICES, STRUCTURE_TYPE, ERR, ERROR,)
Sets the structure (sparsity) of the dynamic equations matrices.
Definition: equations_matrices_routines.f90:4090
iso_varying_string::trim
Definition: iso_varying_string.f90:204
equations_set_constants::equations_set_stokes_equation_three_dim_2
integer(intg), parameter equations_set_stokes_equation_three_dim_2
u=tbd
Definition: equations_set_constants.f90:494
stokes_equations_routines::stokes_equationssetspecificationset
subroutine, public stokes_equationssetspecificationset(equationsSet, specification, err, error,)
Sets the equation specification for a Stokes flow equation of a fluid mechanics equations set...
Definition: Stokes_equations_routines.f90:172
equations_matrices_routines::equations_matrices_linear_storage_type_set
subroutine, public equations_matrices_linear_storage_type_set(EQUATIONS_MATRICES, STORAGE_TYPE, ERR, ERROR,)
Sets the storage type (sparsity) of the linear equations matrices.
Definition: equations_matrices_routines.f90:3898
equations_matrices_routines::equationsmatrices_linearstructuretypeset
subroutine, public equationsmatrices_linearstructuretypeset(EQUATIONS_MATRICES, STRUCTURE_TYPE, ERR, ERROR,)
Sets the structure (sparsity) of the linear equations matrices.
Definition: equations_matrices_routines.f90:4160
equations_mapping_routines::equations_mapping_create_finish
subroutine, public equations_mapping_create_finish(EQUATIONS_MAPPING, ERR, ERROR,)
Finishes the process of creating an equations mapping.
Definition: equations_mapping_routines.f90:982
stokes_equations_routines::stokes_pre_solve_ale_update_parameters
subroutine stokes_pre_solve_ale_update_parameters(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Update mesh parameters for three component Laplace problem.
Definition: Stokes_equations_routines.f90:2949
equations_set_constants::equations_set_stokes_equation_three_dim_3
integer(intg), parameter equations_set_stokes_equation_three_dim_3
u=tbd
Definition: equations_set_constants.f90:495
control_loop_routines::control_loop_get
Returns the specified control loop as indexed by the control loop identifier from the control loop ro...
Definition: control_loop_routines.f90:86
stokes_equations_routines
This module handles all Stokes fluid routines.
Definition: Stokes_equations_routines.f90:46
control_loop_routines::control_loop_type_set
subroutine, public control_loop_type_set(CONTROL_LOOP, LOOP_TYPE, ERR, ERROR,)
Sets/changes the control loop type.
Definition: control_loop_routines.f90:1750
equations_set_constants::equations_set_stokes_equation_three_dim_1
integer(intg), parameter equations_set_stokes_equation_three_dim_1
u=tbd
Definition: equations_set_constants.f90:493
equations_set_constants::equations_set_navier_stokes_equation_three_dim_4
integer(intg), parameter equations_set_navier_stokes_equation_three_dim_4
u=tbd
Definition: equations_set_constants.f90:510
equations_set_constants::equations_set_optimised_stokes_subtype
integer(intg), parameter equations_set_optimised_stokes_subtype
Definition: equations_set_constants.f90:163
problem_constants::problem_transient_stokes_subtype
integer(intg), parameter problem_transient_stokes_subtype
Definition: problem_constants.f90:134
equations_routines::equations_set_equations_get
subroutine, public equations_set_equations_get(EQUATIONS_SET, EQUATIONS, ERR, ERROR,)
Gets the equations for an equations set.
Definition: equations_routines.f90:897
stokes_equations_routines::stokes_finite_element_calculate
subroutine, public stokes_finite_element_calculate(EQUATIONS_SET, ELEMENT_NUMBER, ERR, ERROR,)
Calculates the element stiffness matrices and RHS for a Stokes fluid finite element equations set...
Definition: Stokes_equations_routines.f90:1328
types::equations_set_materials_type
Definition: types.f90:1881
problem_constants::problem_static_stokes_subtype
integer(intg), parameter problem_static_stokes_subtype
Definition: problem_constants.f90:132
constants::partial_derivative_first_derivative_map
integer(intg), dimension(4) partial_derivative_first_derivative_map
PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(nic) gives the partial derivative index for the first derivat...
Definition: constants.f90:254
equations_routines::equations_create_finish
subroutine, public equations_create_finish(EQUATIONS, ERR, ERROR,)
Finish the creation of equations.
Definition: equations_routines.f90:159
domain_mappings
This module handles all domain mappings routines.
Definition: domain_mappings.f90:45
problem_constants::problem_setup_finish_action
integer(intg), parameter problem_setup_finish_action
Finish setup action.
Definition: problem_constants.f90:274
equations_mapping_routines
This module handles all equations mapping routines.
Definition: equations_mapping_routines.f90:45
types::solver_equations_type
Contains information about the solver equations for a solver.
Definition: types.f90:2452
matrix_vector::matrix_compressed_row_storage_type
integer(intg), parameter, public matrix_compressed_row_storage_type
Matrix compressed row storage type.
Definition: matrix_vector.f90:167
equations_matrices_routines::equations_matrices_dynamic_storage_type_set
subroutine, public equations_matrices_dynamic_storage_type_set(EQUATIONS_MATRICES, STORAGE_TYPE, ERR, ERROR,)
Sets the storage type (sparsity) of the dynamic equations matrices.
Definition: equations_matrices_routines.f90:3824
equations_matrices_routines::equations_matrix_diagonal_structure
integer(intg), parameter, public equations_matrix_diagonal_structure
Diagonal matrix structure.
Definition: equations_matrices_routines.f90:73
equations_set_constants::equations_set_gfv_solution_method
integer(intg), parameter equations_set_gfv_solution_method
Grid-based Finite Volume solution method.
Definition: equations_set_constants.f90:354
equations_set_constants::equations_set_setup_geometry_type
integer(intg), parameter equations_set_setup_geometry_type
Geometry setup.
Definition: equations_set_constants.f90:323
constants::global_deriv_s1_s2
integer(intg), parameter global_deriv_s1_s2
Global Cross derivative in the s1 and s2 direction i.e., d^2u/ds1ds2.
Definition: constants.f90:216
types::problem_type
Contains information for a problem.
Definition: types.f90:3221
solver_routines::solver_progress_output
integer(intg), parameter, public solver_progress_output
Progress output from solver routines.
Definition: solver_routines.f90:350
stokes_equations_routines::stokes_post_solve_output_data
subroutine stokes_post_solve_output_data(CONTROL_LOOP, SOLVER, err, error,)
Output data post solve.
Definition: Stokes_equations_routines.f90:3078
equations_set_constants::equations_linear
integer(intg), parameter equations_linear
The equations are linear.
Definition: equations_set_constants.f90:365
types::field_interpolated_point_ptr_type
Definition: types.f90:1136
analytic_analysis_routines::analyticanalysis_output
subroutine, public analyticanalysis_output(FIELD, FILENAME, ERR, ERROR,)
Output the analytic error analysis for a dependent field compared to the analytic values parameter se...
Definition: analytic_analysis_routines.f90:114
types::domain_nodes_type
Contains the topology information for the nodes of a domain.
Definition: types.f90:713
equations_matrices_routines::equations_matrices_create_finish
subroutine, public equations_matrices_create_finish(EQUATIONS_MATRICES, ERR, ERROR,)
Finishes the creation of the equations matrices and RHS for the the equations.
Definition: equations_matrices_routines.f90:404
distributed_matrix_vector
This module handles all distributed matrix vector routines.
Definition: distributed_matrix_vector.f90:45
problem_constants::problem_optimised_stokes_subtype
integer(intg), parameter problem_optimised_stokes_subtype
Definition: problem_constants.f90:135
constants::global_deriv_s1
integer(intg), parameter global_deriv_s1
First global derivative in the s1 direction i.e., du/ds1.
Definition: constants.f90:214
boundary_conditions_routines
This module handles all boundary conditions routines.
Definition: boundary_condition_routines.f90:45
solver_routines
This module handles all solver routines.
Definition: solver_routines.f90:45
equations_mapping_routines::equations_mapping_create_start
subroutine, public equations_mapping_create_start(EQUATIONS, EQUATIONS_MAPPING, ERR, ERROR,)
Finishes the process of creating an equations mapping for a equations set equations.
Definition: equations_mapping_routines.f90:1112
equations_matrices_routines::equations_matrix_unlumped
integer(intg), parameter, public equations_matrix_unlumped
The matrix is not lumped.
Definition: equations_matrices_routines.f90:82
types::equations_matrix_type
Contains information about an equations matrix.
Definition: types.f90:1429
types::quadrature_scheme_type
Contains information for a particular quadrature scheme.
Definition: types.f90:141
field_io_routines
Implements lists of Field IO operation.
Definition: field_IO_routines.f90:46
matrix_vector
This module contains all routines dealing with (non-distributed) matrix and vectors types...
Definition: matrix_vector.f90:129
distributed_matrix_vector::distributed_matrix_block_storage_type
integer(intg), parameter, public distributed_matrix_block_storage_type
Distributed matrix block storage type.
Definition: distributed_matrix_vector.f90:97
equations_matrices_routines::equations_matrix_lumped
integer(intg), parameter, public equations_matrix_lumped
The matrix is "mass" lumped.
Definition: equations_matrices_routines.f90:83
equations_routines::equations_linearity_type_set
subroutine, public equations_linearity_type_set(EQUATIONS, LINEARITY_TYPE, ERR, ERROR,)
Sets/changes the linearity type for equations.
Definition: equations_routines.f90:539
control_loop_routines::control_loop_create_start
subroutine, public control_loop_create_start(PROBLEM, CONTROL_LOOP, ERR, ERROR,)
Start the process of creating a control loop for a problem.
Definition: control_loop_routines.f90:194
equations_set_constants::equations_set_stokes_equation_two_dim_1
integer(intg), parameter equations_set_stokes_equation_two_dim_1
u=tbd
Definition: equations_set_constants.f90:488
problem_constants::problem_setup_solver_equations_type
integer(intg), parameter problem_setup_solver_equations_type
Solver equations setup for a problem.
Definition: problem_constants.f90:265
equations_set_constants::equations_set_pgm_stokes_subtype
integer(intg), parameter equations_set_pgm_stokes_subtype
Definition: equations_set_constants.f90:165
boundary_conditions_routines::boundary_conditions_set_local_dof
Sets a boundary condition on the specified local DOF.
Definition: boundary_condition_routines.f90:143
types::solver_mapping_type
Contains information on the solver mapping between the global equation sets and the solver matrices...
Definition: types.f90:3091
solver_routines::solver_dynamic_scheme_set
subroutine, public solver_dynamic_scheme_set(SOLVER, SCHEME, ERR, ERROR,)
Sets/changes the scheme for a dynamic solver.
Definition: solver_routines.f90:6384
equations_set_constants::equations_set_stokes_equation_two_dim_5
integer(intg), parameter equations_set_stokes_equation_two_dim_5
u=tbd
Definition: equations_set_constants.f90:492
equations_set_constants::equations_set_navier_stokes_equation_three_dim_1
integer(intg), parameter equations_set_navier_stokes_equation_three_dim_1
u=tbd
Definition: equations_set_constants.f90:507
types::field_variable_type
Contains information for a field variable defined on a field.
Definition: types.f90:1289
equations_set_constants::equations_set_fd_solution_method
integer(intg), parameter equations_set_fd_solution_method
Finite Difference solution method.
Definition: equations_set_constants.f90:350
equations_matrices_routines::equations_matrices_sparse_matrices
integer(intg), parameter, public equations_matrices_sparse_matrices
Use sparse equations matrices.
Definition: equations_matrices_routines.f90:90
problem_constants::problem_pgm_stokes_subtype
integer(intg), parameter problem_pgm_stokes_subtype
Definition: problem_constants.f90:137
types::equations_set_setup_type
Contains information on the setup information for an equations set.
Definition: types.f90:1866
types::domain_type
A pointer to the domain decomposition for this domain.
Definition: types.f90:938
problem_constants::problem_setup_start_action
integer(intg), parameter problem_setup_start_action
Start setup action.
Definition: problem_constants.f90:273
solver_routines::solver_equations_time_dependence_type_set
subroutine, public solver_equations_time_dependence_type_set(SOLVER_EQUATIONS, TIME_DEPENDENCE_TYPE, ERR, ERROR,)
Sets/changes the time dependence type for solver equations.
Definition: solver_routines.f90:7788
control_loop_routines
This module handles all control loop routines.
Definition: control_loop_routines.f90:45
solver_routines::solver_cmiss_library
integer(intg), parameter, public solver_cmiss_library
CMISS (internal) solver library.
Definition: solver_routines.f90:103
boundary_conditions_routines::boundary_condition_fixed
integer(intg), parameter, public boundary_condition_fixed
The dof is fixed as a boundary condition.
Definition: boundary_condition_routines.f90:95
base_routines::errors
subroutine, public errors(NAME, ERR, ERROR)
Records the exiting error of the subroutine.
Definition: base_routines.f90:481
equations_set_constants
This module defines all constants shared across equations set routines.
Definition: equations_set_constants.f90:45
equations_set_constants::equations_set_bem_solution_method
integer(intg), parameter equations_set_bem_solution_method
Boundary Element Method solution method.
Definition: equations_set_constants.f90:349
solver_routines::solver_solver_equations_get
subroutine, public solver_solver_equations_get(SOLVER, SOLVER_EQUATIONS, ERR, ERROR,)
Returns a pointer to the solver equations for a solver.
Definition: solver_routines.f90:21230
fluid_mechanics_io_routines::fluid_mechanics_io_read_boundary_conditions
subroutine, public fluid_mechanics_io_read_boundary_conditions(SOLVER_TYPE, BOUNDARY_VALUES, NUMBER_OF_DIMENSIONS, BOUNDARY_CONDITION, OPTION, TIME_STEP, TIME, LENGTH_SCALE)
Reads boundary conditions from a file.
Definition: fluid_mechanics_IO_routines.f90:4128
boundary_conditions_routines::boundary_conditions_variable_get
subroutine, public boundary_conditions_variable_get(BOUNDARY_CONDITIONS, FIELD_VARIABLE, BOUNDARY_CONDITIONS_VARIABLE, ERR, ERROR,)
Find the boundary conditions variable for a given field variable.
Definition: boundary_condition_routines.f90:3635
types::basis_type
Contains all information about a basis .
Definition: types.f90:184
equations_set_constants::equations_set_fv_solution_method
integer(intg), parameter equations_set_fv_solution_method
Finite Volume solution method.
Definition: equations_set_constants.f90:351
matrix_vector::matrix_block_storage_type
integer(intg), parameter, public matrix_block_storage_type
Matrix block storage type.
Definition: matrix_vector.f90:163
solver_routines::solver_dynamic_first_degree
integer(intg), parameter, public solver_dynamic_first_degree
Dynamic solver uses a first degree polynomial for time interpolation.
Definition: solver_routines.f90:271
stokes_equations_routines::stokes_boundaryconditionsanalyticcalculate
subroutine, public stokes_boundaryconditionsanalyticcalculate(EQUATIONS_SET, BOUNDARY_CONDITIONS, ERR, ERROR,)
Calculates the analytic solution and sets the boundary conditions for an analytic problem...
Definition: Stokes_equations_routines.f90:3211
base_routines::flag_error
Flags an error condition.
Definition: base_routines.f90:226
equations_set_constants::equations_set_setup_initial_type
integer(intg), parameter equations_set_setup_initial_type
Initial setup.
Definition: equations_set_constants.f90:322
control_loop_routines::control_loop_create_finish
recursive subroutine, public control_loop_create_finish(CONTROL_LOOP, ERR, ERROR,)
Finish the process of creating a control loop.
Definition: control_loop_routines.f90:153
equations_set_constants::equations_set_setup_analytic_type
integer(intg), parameter equations_set_setup_analytic_type
Analytic setup.
Definition: equations_set_constants.f90:329
base_routines::flagerror
Flags an error condition.
Definition: base_routines.f90:240
solver_routines::solver_linear_type
integer(intg), parameter, public solver_linear_type
A linear solver.
Definition: solver_routines.f90:88
problem_constants::problem_fluid_mechanics_class
integer(intg), parameter problem_fluid_mechanics_class
Definition: problem_constants.f90:56
equations_set_constants::equations_set_stokes_equation_type
integer(intg), parameter equations_set_stokes_equation_type
Definition: equations_set_constants.f90:71
types::equations_matrices_rhs_type
Contains information of the RHS vector for equations matrices.
Definition: types.f90:1500
field_io_routines::field_io_elements_export
subroutine, public field_io_elements_export(FIELDS, FILE_NAME, METHOD, ERR, ERROR,)
Export elemental information into multiple files.
Definition: field_IO_routines.f90:6207
types::field_interpolation_parameters_ptr_type
Definition: types.f90:1151
stokes_equations_routines::stokes_problemspecificationset
subroutine, public stokes_problemspecificationset(problem, problemSpecification, err, error,)
Sets the problem specification for a Stokes fluid problem.
Definition: Stokes_equations_routines.f90:925
equations_set_constants::equations_set_stokes_equation_three_dim_5
integer(intg), parameter equations_set_stokes_equation_three_dim_5
u=tbd
Definition: equations_set_constants.f90:497
distributed_matrix_vector::distributed_matrix_diagonal_storage_type
integer(intg), parameter, public distributed_matrix_diagonal_storage_type
Distributed matrix diagonal storage type.
Definition: distributed_matrix_vector.f90:98
equations_set_constants::equations_set_ale_stokes_subtype
integer(intg), parameter equations_set_ale_stokes_subtype
Definition: equations_set_constants.f90:164
types::equations_mapping_linear_type
Contains information for mapping field variables to the linear matrices in the equations set of the m...
Definition: types.f90:1587
kinds
This module contains all kind definitions.
Definition: kinds.f90:45
fluid_mechanics_io_routines
Temporary IO routines for fluid mechanics.
Definition: fluid_mechanics_IO_routines.f90:47
field_io_routines::field_io_nodes_export
subroutine, public field_io_nodes_export(FIELDS, FILE_NAME, METHOD, ERR, ERROR,)
Export nodal information.
Definition: field_IO_routines.f90:6163
equations_set_constants::equations_set_setup_finish_action
integer(intg), parameter equations_set_setup_finish_action
Finish setup action.
Definition: equations_set_constants.f90:339
distributed_matrix_vector::distributed_matrix_compressed_row_storage_type
integer(intg), parameter, public distributed_matrix_compressed_row_storage_type
Distributed matrix compressed row storage type.
Definition: distributed_matrix_vector.f90:101
types::equations_matrices_dynamic_type
Contains information of the dynamic matrices for equations matrices.
Definition: types.f90:1471
input_output
This module handles all formating and input and output.
Definition: input_output.f90:45