OpenCMISS-Iron Internal API Documentation
Darcy_equations_routines.f90
Go to the documentation of this file.
1 
43 
45 
46 MODULE darcy_equations_routines
47 
48  USE base_routines
49  USE basis_routines
50  USE boundary_conditions_routines
51  USE constants
52  USE control_loop_routines
53  USE comp_environment
54  USE coordinate_routines
55  USE distributed_matrix_vector
56  USE domain_mappings
57  USE equations_routines
58  USE equations_mapping_routines
59  USE equations_matrices_routines
60  USE equations_set_constants
61  USE field_routines
62  USE field_io_routines
63  USE finite_elasticity_routines
64  USE fluid_mechanics_io_routines
65  USE input_output
66  USE iso_varying_string
67  USE kinds
68  USE maths
69  USE matrix_vector
70  USE mesh_routines
71  USE node_routines
72  USE problem_constants
73  USE strings
74  USE solver_routines
75  USE timer
76  USE types
77 
78 #include "macros.h"
79 
80 
81  IMPLICIT NONE
82 
83  PUBLIC darcy_equation_equations_set_setup
84  PUBLIC darcy_equationssetspecificationset
85  PUBLIC darcy_equationssetsolutionmethodset
86  PUBLIC darcy_boundaryconditionsanalyticcalculate
87 
88  PUBLIC darcy_equation_problem_setup
89  PUBLIC darcy_problemspecificationset
90 
91  PUBLIC darcy_equation_finite_element_calculate
92 
93  PUBLIC darcy_equation_pre_solve
94  PUBLIC darcy_equation_post_solve
95  PUBLIC darcy_equation_post_solve_output_data
96 
97  PUBLIC darcy_control_time_loop_pre_loop
98 
99  PUBLIC darcy_presolvestorepreviousiterate
100 
101  PUBLIC darcy_equation_monitor_convergence
102 
103  INTEGER(INTG) :: solver_number_solid,solver_number_mat_properties,solver_number_darcy
104  INTEGER(INTG) :: solver_index_solid,solver_index_mat_properties,solver_index_darcy
105 
106  REAL(DP) :: residual_norm_0
107 
108  LOGICAL :: idebug1, idebug2, idebug3
109 
110 
111 CONTAINS
112 
113  !
114  !================================================================================================================================
115  !
116 
118  SUBROUTINE darcy_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)
119 
120  !Argument variables
121  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
122  INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
123  INTEGER(INTG), INTENT(OUT) :: ERR
124  TYPE(varying_string), INTENT(OUT) :: ERROR
125  !Local Variables
126  TYPE(varying_string) :: LOCAL_ERROR
127 
128  enters("Darcy_EquationsSetSolutionMethodSet",err,error,*999)
129 
130  IF(ASSOCIATED(equations_set)) THEN
131  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
132  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
133  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
134  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
135  & err,error,*999)
136  END IF
137  SELECT CASE(equations_set%SPECIFICATION(3))
138  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
139  & equations_set_transient_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
140  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
141  & equations_set_incompressible_elasticity_driven_darcy_subtype,equations_set_multi_compartment_darcy_subtype, &
142  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
143  SELECT CASE(solution_method)
144  CASE(equations_set_fem_solution_method)
145  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
146  CASE(equations_set_bem_solution_method)
147  CALL flagerror("Not implemented.",err,error,*999)
148  CASE(equations_set_fd_solution_method)
149  CALL flagerror("Not implemented.",err,error,*999)
150  CASE(equations_set_fv_solution_method)
151  CALL flagerror("Not implemented.",err,error,*999)
152  CASE(equations_set_gfem_solution_method)
153  CALL flagerror("Not implemented.",err,error,*999)
154  CASE(equations_set_gfv_solution_method)
155  CALL flagerror("Not implemented.",err,error,*999)
156  CASE DEFAULT
157  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
158  CALL flagerror(local_error,err,error,*999)
159  END SELECT
160  CASE DEFAULT
161  local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
162  & " is not valid for a Darcy equation type of a fluid mechanics equations set class."
163  CALL flagerror(local_error,err,error,*999)
164  END SELECT
165  ELSE
166  CALL flagerror("Equations set is not associated.",err,error,*999)
167  ENDIF
168 
169  exits("Darcy_EquationsSetSolutionMethodSet")
170  RETURN
171 999 errorsexits("Darcy_EquationsSetSolutionMethodSet",err,error)
172  RETURN 1
173  END SUBROUTINE darcy_equationssetsolutionmethodset
174 
175  !
176  !================================================================================================================================
177  !
178 
180  SUBROUTINE darcy_equation_equations_set_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
181 
182  !Argument variables
183  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
184  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
185  INTEGER(INTG), INTENT(OUT) :: ERR
186  TYPE(varying_string), INTENT(OUT) :: ERROR
187  !Local Variables
188  INTEGER(INTG) :: GEOMETRIC_SCALING_TYPE, GEOMETRIC_MESH_COMPONENT,NUMBER_OF_COMPONENTS,NUMBER_OF_DARCY_COMPONENTS
189  INTEGER(INTG) :: imy_matrix,Ncompartments
190  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
191  TYPE(equations_type), POINTER :: EQUATIONS
192  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
193  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
194  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
195  TYPE(equations_set_equations_set_field_type), POINTER :: EQUATIONS_EQUATIONS_SET_FIELD
196  TYPE(field_type), POINTER :: EQUATIONS_SET_FIELD_FIELD
197  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
198  TYPE(equations_set_source_type), POINTER :: EQUATIONS_SOURCE
199  TYPE(varying_string) :: LOCAL_ERROR
200  INTEGER:: DEPENDENT_FIELD_NUMBER_OF_VARIABLES, DEPENDENT_FIELD_NUMBER_OF_COMPONENTS
201  INTEGER:: DEPENDENT_FIELD_ELASTICITY_NUMBER_OF_COMPONENTS, DEPENDENT_FIELD_DARCY_NUMBER_OF_COMPONENTS
202  INTEGER:: INDEPENDENT_FIELD_NUMBER_OF_VARIABLES, INDEPENDENT_FIELD_NUMBER_OF_COMPONENTS
203  INTEGER:: NUMBER_OF_DIMENSIONS, GEOMETRIC_COMPONENT_NUMBER
204  INTEGER:: MATERIAL_FIELD_NUMBER_OF_VARIABLES, MATERIAL_FIELD_NUMBER_OF_COMPONENTS
205  INTEGER:: MESH_COMPONENT,MATERIAL_FIELD_NUMBER_OF_U_VAR_COMPONENTS,MATERIAL_FIELD_NUMBER_OF_V_VAR_COMPONENTS, &
206  & MATERIAL_FIELD_NUMBER_OF_U1_VAR_COMPONENTS
207  INTEGER:: i,component_idx
208 
209  INTEGER(INTG) :: num_var,num_var_count
210  INTEGER(INTG) :: EQUATIONS_SET_FIELD_NUMBER_OF_VARIABLES,EQUATIONS_SET_FIELD_NUMBER_OF_COMPONENTS,NUMBER_OF_SOURCE_COMPONENTS
211  INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
212  INTEGER(INTG), ALLOCATABLE :: VARIABLE_TYPES(:),VARIABLE_U_TYPES(:),COUPLING_MATRIX_STORAGE_TYPE(:), &
213  & COUPLING_MATRIX_STRUCTURE_TYPE(:)
214 
215  enters("DARCY_EQUATION_EQUATIONS_SET_SETUP",err,error,*999)
216 
217  NULLIFY(equations)
218  NULLIFY(equations_mapping)
219  NULLIFY(equations_matrices)
220  NULLIFY(geometric_decomposition)
221  NULLIFY(equations_equations_set_field)
222  NULLIFY(equations_set_field_field)
223  NULLIFY(equations_set_field_data)
224 
225  IF(ASSOCIATED(equations_set)) THEN
226  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
227  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
228  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
229  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
230  & err,error,*999)
231  END IF
232  SELECT CASE(equations_set%SPECIFICATION(3))
233  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
234  & equations_set_transient_darcy_subtype, equations_set_incompressible_finite_elasticity_darcy_subtype, &
235  & equations_set_transient_ale_darcy_subtype, equations_set_elasticity_darcy_inria_model_subtype, &
236  & equations_set_incompressible_elasticity_driven_darcy_subtype,equations_set_multi_compartment_darcy_subtype, &
237  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
238  SELECT CASE(equations_set_setup%SETUP_TYPE)
239 
240  !-----------------------------------------------------------------
241  ! s o l u t i o n m e t h o d
242  !-----------------------------------------------------------------
243  CASE(equations_set_setup_initial_type)
244  SELECT CASE(equations_set%SPECIFICATION(3))
245  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
246  & equations_set_transient_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
247  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
248  & equations_set_incompressible_elasticity_driven_darcy_subtype)
249  SELECT CASE(equations_set_setup%ACTION_TYPE)
250  CASE(equations_set_setup_start_action)
251  CALL darcy_equationssetsolutionmethodset(equations_set,equations_set_fem_solution_method,err,error,*999)
252  CASE(equations_set_setup_finish_action)
253  !do nothing
254  CASE DEFAULT
255  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
256  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
257  & " is invalid for a standard or quasistatic Darcy equation."
258  CALL flagerror(local_error,err,error,*999)
259  END SELECT
260  CASE(equations_set_multi_compartment_darcy_subtype,equations_set_incompressible_elast_multi_comp_darcy_subtype)
261  SELECT CASE(equations_set_setup%ACTION_TYPE)
262  CASE(equations_set_setup_start_action)
263  CALL darcy_equationssetsolutionmethodset(equations_set,equations_set_fem_solution_method,err,error,*999)
264 
265  equations_set_field_number_of_variables = 1
266  equations_set_field_number_of_components = 2
267  equations_equations_set_field=>equations_set%EQUATIONS_SET_FIELD
268  IF(equations_equations_set_field%EQUATIONS_SET_FIELD_AUTO_CREATED) THEN
269  !Create the auto created equations set field
270  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
271  & equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,err,error,*999)
272  equations_set_field_field=>equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD
273  CALL field_label_set(equations_set_field_field,"Equations Set Field",err,error,*999)
274  CALL field_type_set_and_lock(equations_set_field_field,field_general_type,&
275  & err,error,*999)
276  CALL field_dependent_type_set_and_lock(equations_set_field_field,&
277  & field_independent_type,err,error,*999)
278  CALL field_number_of_variables_set(equations_set_field_field, &
279  & equations_set_field_number_of_variables,err,error,*999)
280  CALL field_variable_types_set_and_lock(equations_set_field_field,&
281  & [field_u_variable_type],err,error,*999)
282  CALL field_dimension_set_and_lock(equations_set_field_field,field_u_variable_type, &
283  & field_vector_dimension_type,err,error,*999)
284  CALL field_data_type_set_and_lock(equations_set_field_field,field_u_variable_type, &
285  & field_intg_type,err,error,*999)
286  CALL field_number_of_components_set_and_lock(equations_set_field_field,&
287  & field_u_variable_type,equations_set_field_number_of_components,err,error,*999)
288  ELSE
289  !Check the user specified field
290  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
291  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
292  CALL field_number_of_variables_check(equations_set_setup%FIELD,equations_set_field_number_of_variables, &
293  & err,error,*999)
294  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
295  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
296  & err,error,*999)
297  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_intg_type,err,error,*999)
298  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
299  & equations_set_field_number_of_components,err,error,*999)
300  ENDIF
301  CASE(equations_set_setup_finish_action)
302  IF(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_AUTO_CREATED) THEN
303  CALL field_create_finish(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,err,error,*999)
304  CALL field_component_values_initialise(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,&
305  & field_u_variable_type,field_values_set_type, 1, 1_intg, err, error, *999)
306  CALL field_component_values_initialise(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,&
307  & field_u_variable_type,field_values_set_type, 2, 1_intg, err, error, *999)
308  ENDIF
309 !!TODO: Check valid setup
310  CASE DEFAULT
311  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
312  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
313  & " is invalid for a standard or quasistatic Darcy equation."
314  CALL flagerror(local_error,err,error,*999)
315  END SELECT
316  END SELECT
317 
318  !-----------------------------------------------------------------
319  ! g e o m e t r y f i e l d
320  !-----------------------------------------------------------------
321  CASE(equations_set_setup_geometry_type)
322  SELECT CASE(equations_set%SPECIFICATION(3))
323  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, &
324  & equations_set_transient_darcy_subtype)
325  !Do nothing
326  CASE(equations_set_ale_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
327  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
328  & equations_set_incompressible_elasticity_driven_darcy_subtype,equations_set_multi_compartment_darcy_subtype, &
329  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
330 
331  SELECT CASE(equations_set_setup%ACTION_TYPE)
332  CASE(equations_set_setup_start_action)
333 
334  field_variable=>equations_set%GEOMETRY%GEOMETRIC_FIELD%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
335 
336  CALL field_parametersetensurecreated(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
337  & field_initial_values_set_type, err, error, *999)
338 
339  CALL field_parametersetensurecreated(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
340  & field_previous_values_set_type, err, error, *999)
341 
342  CALL field_parametersetensurecreated(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
343  & field_mesh_displacement_set_type, err, error, *999)
344 
345  CALL field_parametersetensurecreated(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
346  & field_mesh_velocity_set_type, err, error, *999)
347 
348  CALL field_parametersetensurecreated(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
349  & field_negative_mesh_velocity_set_type, err, error, *999)
350 
351  IF(equations_set%SPECIFICATION(3)==equations_set_multi_compartment_darcy_subtype .OR. &
352  equations_set%SPECIFICATION(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
353  !Create the equations set field for multi-compartment Darcy
354  equations_set_field_number_of_components = 2
355 
356  equations_equations_set_field=>equations_set%EQUATIONS_SET_FIELD
357  equations_set_field_field=>equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD
358 
359  IF(equations_equations_set_field%EQUATIONS_SET_FIELD_AUTO_CREATED) THEN
360  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
361  CALL field_mesh_decomposition_set_and_lock(equations_set_field_field,&
362  & geometric_decomposition,err,error,*999)
363  CALL field_geometric_field_set_and_lock(equations_set_field_field,&
364  & equations_set%GEOMETRY%GEOMETRIC_FIELD,err,error,*999)
365  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
366  & 1,geometric_component_number,err,error,*999)
367  DO component_idx = 1, equations_set_field_number_of_components
368  CALL field_component_mesh_component_set_and_lock(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD, &
369  & field_u_variable_type,component_idx,geometric_component_number,err,error,*999)
370  CALL field_component_interpolation_set_and_lock(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD, &
371  & field_u_variable_type,component_idx,field_constant_interpolation,err,error,*999)
372  END DO
373 
374  !Default the field scaling to that of the geometric field
375  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
376  CALL field_scaling_type_set(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,geometric_scaling_type, &
377  & err,error,*999)
378  ELSE
379  !Do nothing
380  ENDIF
381  ENDIF
382  CASE(equations_set_setup_finish_action)
383  ! do nothing
384  CASE DEFAULT
385  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
386  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
387  & " is invalid for a linear diffusion equation."
388  CALL flagerror(local_error,err,error,*999)
389  END SELECT
390  END SELECT
391 
392  !-----------------------------------------------------------------
393  ! d e p e n d e n t f i e l d
394  !-----------------------------------------------------------------
395  CASE(equations_set_setup_dependent_type)
396  SELECT CASE(equations_set%SPECIFICATION(3))
397  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
398  & equations_set_transient_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
399  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
400  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
401  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
402  SELECT CASE(equations_set_setup%ACTION_TYPE)
403  CASE(equations_set_setup_start_action)
404  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
405  !Create the auto created dependent field
406  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
407  & equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
408  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_geometric_general_type,err,error,*999)
409  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
410 
411  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
412  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
413  & err,error,*999)
414  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
415  & equations_set%GEOMETRY%GEOMETRIC_FIELD,err,error,*999)
416 
417  dependent_field_number_of_variables = 2 ! U and the normal component of its flux
418  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
419  & dependent_field_number_of_variables,err,error,*999)
420  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,[field_u_variable_type, &
421  & field_deludeln_variable_type],err,error,*999)
422  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,"U",err,error,*999)
423  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,"del U/del n", &
424  & err,error,*999)
425  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
426  & field_vector_dimension_type,err,error,*999)
427  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
428  & field_vector_dimension_type,err,error,*999)
429  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
430  & field_dp_type,err,error,*999)
431  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
432  & field_dp_type,err,error,*999)
433 
434  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
435  & number_of_dimensions, err, error, *999)
436  dependent_field_number_of_components = number_of_dimensions + 1
437  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, field_u_variable_type, &
438  & dependent_field_number_of_components, err, error, *999)
439  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
440  & field_deludeln_variable_type, dependent_field_number_of_components, err, error, *999)
441  !Default to the geometric interpolation setup
442  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
443  & geometric_mesh_component,err,error,*999)
444  DO i=1,dependent_field_number_of_components
445  IF( i < dependent_field_number_of_components ) THEN
446  !Set velocity mesh component (default to the geometric one)
447  mesh_component = geometric_mesh_component
448  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
449  & i, mesh_component,err,error,*999)
450  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
451  & i, mesh_component,err,error,*999)
452  ELSE
453  !Set pressure mesh component (default to the geometric one)
454  mesh_component = geometric_mesh_component
455  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
456  & i, mesh_component,err,error,*999)
457  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
458  & i, mesh_component,err,error,*999)
459  ENDIF
460  ENDDO
461 
462  SELECT CASE(equations_set%SOLUTION_METHOD)
463  CASE(equations_set_fem_solution_method)
464  DO i = 1, dependent_field_number_of_components
465  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
466  & field_u_variable_type,i,field_node_based_interpolation,err,error,*999)
467  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
468  & field_deludeln_variable_type,i,field_node_based_interpolation,err,error,*999)
469  ENDDO
470  !Default the scaling to the geometric field scaling
471  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
472  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
473  CASE(equations_set_bem_solution_method)
474  CALL flagerror("Not implemented.",err,error,*999)
475  CASE(equations_set_fd_solution_method)
476  CALL flagerror("Not implemented.",err,error,*999)
477  CASE(equations_set_fv_solution_method)
478  CALL flagerror("Not implemented.",err,error,*999)
479  CASE(equations_set_gfem_solution_method)
480  CALL flagerror("Not implemented.",err,error,*999)
481  CASE(equations_set_gfv_solution_method)
482  CALL flagerror("Not implemented.",err,error,*999)
483  CASE DEFAULT
484  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
485  & " is invalid."
486  CALL flagerror(local_error,err,error,*999)
487  END SELECT
488  ELSE
489  !-----------------------------------
490  ! DEPENDENT_FIELD: not AUTO_CREATED
491  !-----------------------------------
492  SELECT CASE(equations_set%SPECIFICATION(3))
493  CASE(equations_set_incompressible_finite_elasticity_darcy_subtype, &
494  & equations_set_elasticity_darcy_inria_model_subtype, &
495  & equations_set_incompressible_elasticity_driven_darcy_subtype)
496  !-----------------------------------------------------------------------
497  ! Check the shared dependent field set up in finite elasticity routines
498  !-----------------------------------------------------------------------
499  CALL field_type_check(equations_set_setup%FIELD,field_geometric_general_type,err,error,*999)
500  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
501  CALL field_number_of_variables_check(equations_set_setup%FIELD,4,err,error,*999)
502  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type,field_deludeln_variable_type, &
503  & field_v_variable_type,field_delvdeln_variable_type],err,error,*999)
504  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
505  & err,error,*999)
506  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_vector_dimension_type, &
507  & err,error,*999)
508  CALL field_dimension_check(equations_set_setup%FIELD,field_v_variable_type,field_vector_dimension_type, &
509  & err,error,*999)
510  CALL field_dimension_check(equations_set_setup%FIELD,field_delvdeln_variable_type,field_vector_dimension_type, &
511  & err,error,*999)
512  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
513  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
514  CALL field_data_type_check(equations_set_setup%FIELD,field_v_variable_type,field_dp_type,err,error,*999)
515  CALL field_data_type_check(equations_set_setup%FIELD,field_delvdeln_variable_type,field_dp_type,err,error,*999)
516  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
517  & number_of_dimensions, err, error, *999)
518 
519  SELECT CASE(equations_set%SPECIFICATION(3))
520  CASE(equations_set_elasticity_darcy_inria_model_subtype) !compressible elasticity
521  dependent_field_elasticity_number_of_components = number_of_dimensions
522  dependent_field_darcy_number_of_components = number_of_dimensions + 2 !(u,v,w,p,m)
523  CASE(equations_set_incompressible_finite_elasticity_darcy_subtype)
524  dependent_field_elasticity_number_of_components = number_of_dimensions + 1
525  dependent_field_darcy_number_of_components = number_of_dimensions + 1
526  CASE(equations_set_incompressible_elasticity_driven_darcy_subtype)
527  dependent_field_elasticity_number_of_components = number_of_dimensions + 1 !(u1,u2,u3,p)
528  dependent_field_darcy_number_of_components = number_of_dimensions + 1 !(u,v,w,m)
529  END SELECT
530 
531  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
532  & dependent_field_elasticity_number_of_components,err,error,*999)
533  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
534  & dependent_field_elasticity_number_of_components,err,error,*999)
535  CALL field_number_of_components_check(equations_set_setup%FIELD,field_v_variable_type, &
536  & dependent_field_darcy_number_of_components,err,error,*999)
537  CALL field_number_of_components_check(equations_set_setup%FIELD,field_delvdeln_variable_type, &
538  & dependent_field_darcy_number_of_components,err,error,*999)
539 
540  SELECT CASE(equations_set%SOLUTION_METHOD)
541  CASE(equations_set_fem_solution_method)
542  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
543  & field_node_based_interpolation,err,error,*999)
544  !Mind that elastic hydrostatic pressure might be interpolated element-wise
545  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
546  & field_node_based_interpolation,err,error,*999)
547  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_v_variable_type,1, &
548  & field_node_based_interpolation,err,error,*999)
549  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_delvdeln_variable_type,1, &
550  & field_node_based_interpolation,err,error,*999)
551  CASE(equations_set_bem_solution_method)
552  CALL flagerror("Not implemented.",err,error,*999)
553  CASE(equations_set_fd_solution_method)
554  CALL flagerror("Not implemented.",err,error,*999)
555  CASE(equations_set_fv_solution_method)
556  CALL flagerror("Not implemented.",err,error,*999)
557  CASE(equations_set_gfem_solution_method)
558  CALL flagerror("Not implemented.",err,error,*999)
559  CASE(equations_set_gfv_solution_method)
560  CALL flagerror("Not implemented.",err,error,*999)
561  CASE DEFAULT
562  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
563  & " is invalid."
564  CALL flagerror(local_error,err,error,*999)
565  END SELECT
566  CASE(equations_set_incompressible_elast_multi_comp_darcy_subtype)
567  !-----------------------------------------------------------------------
568  ! Check the shared dependent field set up in finite elasticity routines
569  ! Must have 2+2*Ncompartments number of variable types
570  !-----------------------------------------------------------------------
571  CALL field_type_check(equations_set_setup%FIELD,field_geometric_general_type,err,error,*999)
572  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
573  !Get the number of Darcy compartments from the equations set field
574  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
575  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
576  & field_values_set_type,equations_set_field_data,err,error,*999)
577  ncompartments=equations_set_field_data(2)
578  CALL field_number_of_variables_check(equations_set_setup%FIELD,(2+2*ncompartments),err,error,*999)
579  ALLOCATE(variable_types(2*ncompartments+2))
580  DO num_var=1,ncompartments+1
581  variable_types(2*num_var-1)=field_u_variable_type+(field_number_of_variable_subtypes*(num_var-1))
582  variable_types(2*num_var)=field_deludeln_variable_type+(field_number_of_variable_subtypes*(num_var-1))
583  ENDDO
584  CALL field_variable_types_check(equations_set_setup%FIELD,variable_types,err,error,*999)
585 
586  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
587  & number_of_dimensions,err,error,*999)
588  number_of_components=number_of_dimensions+1
589  number_of_darcy_components=number_of_dimensions+1
590 
591  DO num_var=1,2*ncompartments+2
592  CALL field_dimension_check(equations_set_setup%FIELD,variable_types(num_var),field_vector_dimension_type, &
593  & err,error,*999)
594  CALL field_data_type_check(equations_set_setup%FIELD,variable_types(num_var),field_dp_type,err,error,*999)
595  CALL field_number_of_components_check(equations_set_setup%FIELD,variable_types(num_var),number_of_components, &
596  & err,error,*999)
597  ENDDO
598 
599  SELECT CASE(equations_set%SOLUTION_METHOD)
600  CASE(equations_set_fem_solution_method)
601  !Elasticity:
602  DO component_idx=1,number_of_dimensions
603  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,component_idx, &
604  & field_node_based_interpolation,err,error,*999)
605  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,component_idx,&
606  & field_node_based_interpolation,err,error,*999)
607  ENDDO !component_idx
608  !If solid hydrostatic pressure is driving Darcy flow, check that pressure uses node based interpolation
609  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,number_of_components,&
610  & field_node_based_interpolation,err,error,*999)
611  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
612  & number_of_components,field_node_based_interpolation,err,error,*999)
613  DO num_var=3,2*ncompartments+2
614  !Darcy:
615  DO component_idx=1,number_of_darcy_components
616  CALL field_component_interpolation_check(equations_set_setup%FIELD,variable_types(num_var),component_idx, &
617  & field_node_based_interpolation,err,error,*999)
618  ENDDO !component_idx
619  ENDDO
620  CASE(equations_set_bem_solution_method)
621  CALL flagerror("Not implemented.",err,error,*999)
622  CASE(equations_set_fd_solution_method)
623  CALL flagerror("Not implemented.",err,error,*999)
624  CASE(equations_set_fv_solution_method)
625  CALL flagerror("Not implemented.",err,error,*999)
626  CASE(equations_set_gfem_solution_method)
627  CALL flagerror("Not implemented.",err,error,*999)
628  CASE(equations_set_gfv_solution_method)
629  CALL flagerror("Not implemented.",err,error,*999)
630  CASE DEFAULT
631  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
632  & " is invalid."
633  CALL flagerror(local_error,err,error,*999)
634  END SELECT
635  CASE(equations_set_multi_compartment_darcy_subtype)
636  !Check the field created by Darcy routines for the multi-compartment model
637  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
638  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
639  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
640  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
641  & field_values_set_type,equations_set_field_data,err,error,*999)
642  ncompartments=equations_set_field_data(2)
643  CALL field_number_of_variables_check(equations_set_setup%FIELD,2*ncompartments,err,error,*999)
644  !Create & populate array storing all of the relevant variable types against which to check the field variables
645  ALLOCATE(variable_types(2*ncompartments))
646  DO num_var=1,ncompartments
647  variable_types(2*num_var-1)=field_u_variable_type+(field_number_of_variable_subtypes*(num_var-1))
648  variable_types(2*num_var)=field_deludeln_variable_type+(field_number_of_variable_subtypes*(num_var-1))
649  ENDDO
650  CALL field_variable_types_check(equations_set_setup%FIELD,variable_types,err,error,*999)
651 
652  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
653  & number_of_dimensions,err,error,*999)
654  DO num_var=1,2*ncompartments
655  CALL field_dimension_check(equations_set_setup%FIELD,variable_types(num_var), &
656  & field_vector_dimension_type,err,error,*999)
657  CALL field_data_type_check(equations_set_setup%FIELD,variable_types(num_var),field_dp_type,err,error,*999)
658  CALL field_number_of_components_check(equations_set_setup%FIELD,variable_types(num_var), &
659  & number_of_dimensions+1,err,error,*999)
660  ENDDO
661 
662  SELECT CASE(equations_set%SOLUTION_METHOD)
663  CASE(equations_set_fem_solution_method)
664  component_idx=1
665  DO num_var=1,2*ncompartments
666  DO component_idx=1,number_of_dimensions+1
667  CALL field_component_interpolation_check(equations_set_setup%FIELD,variable_types(num_var),component_idx, &
668  & field_node_based_interpolation,err,error,*999)
669  !NOTE-pressure might use element based interpolation - need to account for this
670  ENDDO
671  ENDDO
672  CASE(equations_set_bem_solution_method)
673  CALL flagerror("Not implemented.",err,error,*999)
674  CASE(equations_set_fd_solution_method)
675  CALL flagerror("Not implemented.",err,error,*999)
676  CASE(equations_set_fv_solution_method)
677  CALL flagerror("Not implemented.",err,error,*999)
678  CASE(equations_set_gfem_solution_method)
679  CALL flagerror("Not implemented.",err,error,*999)
680  CASE(equations_set_gfv_solution_method)
681  CALL flagerror("Not implemented.",err,error,*999)
682  CASE DEFAULT
683  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
684  & " is invalid."
685  CALL flagerror(local_error,err,error,*999)
686  END SELECT
687 
688  CASE DEFAULT
689  !--------------------------------
690  ! Check the user specified field
691  !--------------------------------
692  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
693  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
694  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
695  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type,field_deludeln_variable_type],&
696  & err,error,*999)
697  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
698  & err,error,*999)
699  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_vector_dimension_type, &
700  & err,error,*999)
701  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
702  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
703  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
704  & number_of_dimensions, err, error, *999)
705  dependent_field_number_of_components = number_of_dimensions + 1
706  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
707  & dependent_field_number_of_components,err,error,*999)
708  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
709  & dependent_field_number_of_components,err,error,*999)
710 
711  SELECT CASE(equations_set%SOLUTION_METHOD)
712  CASE(equations_set_fem_solution_method)
713  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
714  & field_node_based_interpolation,err,error,*999)
715  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
716  & field_node_based_interpolation,err,error,*999)
717  CASE(equations_set_bem_solution_method)
718  CALL flagerror("Not implemented.",err,error,*999)
719  CASE(equations_set_fd_solution_method)
720  CALL flagerror("Not implemented.",err,error,*999)
721  CASE(equations_set_fv_solution_method)
722  CALL flagerror("Not implemented.",err,error,*999)
723  CASE(equations_set_gfem_solution_method)
724  CALL flagerror("Not implemented.",err,error,*999)
725  CASE(equations_set_gfv_solution_method)
726  CALL flagerror("Not implemented.",err,error,*999)
727  CASE DEFAULT
728  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
729  & " is invalid."
730  CALL flagerror(local_error,err,error,*999)
731  END SELECT
732  END SELECT ! on (EQUATIONS_SET%SPECIFICATION(3))
733  ENDIF ! on (EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED)
734  CASE(equations_set_setup_finish_action)
735  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
736  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
737  CALL field_parameter_set_create(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
738  & field_initial_values_set_type,err,error,*999)
739  ENDIF
740  IF(equations_set%SPECIFICATION(3)/=equations_set_incompressible_elast_multi_comp_darcy_subtype)THEN
741  !Actually, only needed for PGM (for elasticity_Darcy defined in elasticity V var):
742  CALL field_parameter_set_create(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
743  & field_relative_velocity_set_type,err,error,*999)
744  ENDIF
745  CASE DEFAULT
746  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
747  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
748  & " is invalid for a standard, quasistatic or ALE Darcy equation"
749  CALL flagerror(local_error,err,error,*999)
750  END SELECT
751  END SELECT
752 
753  !-----------------------------------------------------------------
754  ! I N d e p e n d e n t f i e l d
755  !-----------------------------------------------------------------
756  CASE(equations_set_setup_independent_type)
757  SELECT CASE(equations_set%SPECIFICATION(3))
758  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
759  & equations_set_transient_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
760  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
761  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
762  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
763  !\todo: revise: do they all need an independent field ?
764  SELECT CASE(equations_set_setup%ACTION_TYPE)
765  CASE(equations_set_setup_start_action)
766  IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
767  !Create the auto created INdependent field
768  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
769  & equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
770  CALL field_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_general_type,err,error,*999)
771 
772  CALL field_dependent_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_independent_type,&
773  & err,error,*999)
774 
775  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
776  CALL field_mesh_decomposition_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,geometric_decomposition, &
777  & err,error,*999)
778  CALL field_geometric_field_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
779  & equations_set%GEOMETRY%GEOMETRIC_FIELD,err,error,*999)
780 
781  independent_field_number_of_variables = 2 ! U and the normal component of its flux
782  CALL field_number_of_variables_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
783  & independent_field_number_of_variables,err,error,*999)
784  CALL field_variable_types_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,[field_u_variable_type, &
785  & field_deludeln_variable_type],err,error,*999)
786  CALL field_variable_label_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,"Independent U", &
787  & err,error,*999)
788  CALL field_variable_label_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_deludeln_variable_type, &
789  & "Independent del U/del n",err,error,*999)
790  CALL field_dimension_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
791  & field_vector_dimension_type,err,error,*999)
792  CALL field_dimension_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_deludeln_variable_type, &
793  & field_vector_dimension_type,err,error,*999)
794  CALL field_data_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
795  & field_dp_type,err,error,*999)
796  CALL field_data_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_deludeln_variable_type, &
797  & field_dp_type,err,error,*999)
798 
799  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
800  & number_of_dimensions, err, error, *999)
801  independent_field_number_of_components = number_of_dimensions !+ 1
802  CALL field_number_of_components_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, field_u_variable_type, &
803  & independent_field_number_of_components, err, error, *999)
804  CALL field_number_of_components_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
805  & field_deludeln_variable_type, independent_field_number_of_components, err, error, *999)
806  !Default to the geometric interpolation setup
807  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
808  & geometric_mesh_component,err,error,*999)
809  DO i=1,independent_field_number_of_components
810  IF( i < independent_field_number_of_components ) THEN
811  !Set velocity mesh component (default to the geometric one)
812  mesh_component = geometric_mesh_component
813  CALL field_component_mesh_component_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
814  & i, mesh_component,err,error,*999)
815  CALL field_component_mesh_component_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,&
816  & field_deludeln_variable_type, i, mesh_component,err,error,*999)
817  ELSE
818  !Set pressure mesh component (default to the geometric one)
819  mesh_component = geometric_mesh_component
820  CALL field_component_mesh_component_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
821  & i, mesh_component,err,error,*999)
822  CALL field_component_mesh_component_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,&
823  & field_deludeln_variable_type, i, mesh_component,err,error,*999)
824  ENDIF
825  ENDDO
826 
827  SELECT CASE(equations_set%SOLUTION_METHOD)
828  CASE(equations_set_fem_solution_method)
829  DO i = 1, independent_field_number_of_components
830  CALL field_component_interpolation_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
831  & field_u_variable_type,i,field_node_based_interpolation,err,error,*999)
832  CALL field_component_interpolation_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
833  & field_deludeln_variable_type,i,field_node_based_interpolation,err,error,*999)
834  ENDDO
835  !Default the scaling to the geometric field scaling
836  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
837  CALL field_scaling_type_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
838  CASE(equations_set_bem_solution_method)
839  CALL flagerror("Not implemented.",err,error,*999)
840  CASE(equations_set_fd_solution_method)
841  CALL flagerror("Not implemented.",err,error,*999)
842  CASE(equations_set_fv_solution_method)
843  CALL flagerror("Not implemented.",err,error,*999)
844  CASE(equations_set_gfem_solution_method)
845  CALL flagerror("Not implemented.",err,error,*999)
846  CASE(equations_set_gfv_solution_method)
847  CALL flagerror("Not implemented.",err,error,*999)
848  CASE DEFAULT
849  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
850  & " is invalid."
851  CALL flagerror(local_error,err,error,*999)
852  END SELECT
853  ELSE
854  !Check the user specified field
855  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
856 
857  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
858 
859  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
860  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type,field_deludeln_variable_type], &
861  & err,error,*999)
862  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
863  & err,error,*999)
864  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_vector_dimension_type, &
865  & err,error,*999)
866  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
867  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
868  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
869  & number_of_dimensions, err, error, *999)
870  independent_field_number_of_components = number_of_dimensions !+ 1
871  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
872  & independent_field_number_of_components,err,error,*999)
873  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
874  & independent_field_number_of_components,err,error,*999)
875 
876  SELECT CASE(equations_set%SOLUTION_METHOD)
877  CASE(equations_set_fem_solution_method)
878  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
879  & field_node_based_interpolation,err,error,*999)
880  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
881  & field_node_based_interpolation,err,error,*999)
882  CASE(equations_set_bem_solution_method)
883  CALL flagerror("Not implemented.",err,error,*999)
884  CASE(equations_set_fd_solution_method)
885  CALL flagerror("Not implemented.",err,error,*999)
886  CASE(equations_set_fv_solution_method)
887  CALL flagerror("Not implemented.",err,error,*999)
888  CASE(equations_set_gfem_solution_method)
889  CALL flagerror("Not implemented.",err,error,*999)
890  CASE(equations_set_gfv_solution_method)
891  CALL flagerror("Not implemented.",err,error,*999)
892  CASE DEFAULT
893  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
894  & " is invalid."
895  CALL flagerror(local_error,err,error,*999)
896  END SELECT
897  ENDIF
898  CASE(equations_set_setup_finish_action)
899  IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
900  CALL field_create_finish(equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
901  ENDIF
902  CASE DEFAULT
903  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
904  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
905  & " is invalid for a standard, quasistatic or ALE Darcy equation"
906  CALL flagerror(local_error,err,error,*999)
907  END SELECT
908  END SELECT
909 
910  !-----------------------------------------------------------------
911  ! m a t e r i a l f i e l d
912  !-----------------------------------------------------------------
913  CASE(equations_set_setup_materials_type)
914  SELECT CASE(equations_set%SPECIFICATION(3))
915  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
916  & equations_set_transient_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
917  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
918  & equations_set_incompressible_elasticity_driven_darcy_subtype, equations_set_multi_compartment_darcy_subtype)
919  material_field_number_of_variables = 1
920  SELECT CASE(equations_set%SPECIFICATION(3))
921  CASE(equations_set_standard_darcy_subtype,equations_set_quasistatic_darcy_subtype,equations_set_ale_darcy_subtype)
922  !Porosity + scalar permeability/viscosity
923  material_field_number_of_components = 2
924  CASE DEFAULT
925  !Porosity + symmetric permeability/viscosity tensor
926  material_field_number_of_components = 7
927  END SELECT
928  SELECT CASE(equations_set_setup%ACTION_TYPE)
929  CASE(equations_set_setup_start_action)
930  equations_materials=>equations_set%MATERIALS
931  IF(ASSOCIATED(equations_materials)) THEN
932  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
933  !Create the auto created materials field
934  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
935  & materials_field,err,error,*999)
936  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
937  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
938  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
939  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
940  & err,error,*999)
941  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
942  & geometric_field,err,error,*999)
943  CALL field_number_of_variables_set(equations_materials%MATERIALS_FIELD, &
944  & material_field_number_of_variables,err,error,*999)
945  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type], &
946  & err,error,*999)
947  CALL field_variable_label_set(equations_materials%MATERIALS_FIELD,field_u_variable_type,"Material", &
948  & err,error,*999)
949  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
950  & field_vector_dimension_type,err,error,*999)
951  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
952  & field_dp_type,err,error,*999)
953  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
954  & material_field_number_of_components,err,error,*999)
955  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
956  & 1,geometric_component_number,err,error,*999)
957 
958  !Auto-created / default is node_based_interpolation: that's an expensive default ...
959  !Maybe default should be constant; node_based should be requested by the user \todo
960  DO i = 1, material_field_number_of_components
961  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
962  & i,field_node_based_interpolation,err,error,*999)
963  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
964  & i,geometric_component_number,err,error,*999)
965  END DO
966 
967  !Default the field scaling to that of the geometric field
968  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
969  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
970  ELSE
971  !Check the user specified field
972  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
973  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
974  CALL field_number_of_variables_check(equations_set_setup%FIELD,material_field_number_of_variables,err,error,*999)
975  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
976  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
977  & err,error,*999)
978  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
979  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
980  & material_field_number_of_components,err,error,*999)
981  ENDIF
982  ELSE
983  CALL flagerror("Equations set materials is not associated.",err,error,*999)
984  ENDIF
985  CASE(equations_set_setup_finish_action)
986  equations_materials=>equations_set%MATERIALS
987  IF( ASSOCIATED(equations_materials) ) THEN
988  IF( equations_materials%MATERIALS_FIELD_AUTO_CREATED ) THEN
989  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
990  !Set the default values for the materials field
991  DO i=1,material_field_number_of_components
992  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
993  & field_values_set_type, i, 1.0_dp, err, error, *999)
994  ENDDO
995  ENDIF
996  ELSE
997  CALL flagerror("Equations set materials is not associated.",err,error,*999)
998  ENDIF
999  CASE DEFAULT
1000  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1001  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1002  & " is invalid for a standard, quasistatic or ALE Darcy equation."
1003  CALL flagerror(local_error,err,error,*999)
1004  END SELECT
1005  CASE(equations_set_incompressible_elast_multi_comp_darcy_subtype)
1006  !Materials field needs two extra variable types
1007  !The V variable type stores the Darcy coupling coefficients that govern flux between compartments
1008  !The U1 variable type stores the parameters for the constitutive laws that determine the partial pressure in each compartment
1009  !For a first attempt at this, it will be assumed that the functional form of this law is the same for each compartment, with only the paramenters varying (default will be three components)
1010  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1011  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1012  & field_values_set_type,equations_set_field_data,err,error,*999)
1013  ncompartments=equations_set_field_data(2)
1014  material_field_number_of_variables = 3
1015  material_field_number_of_u_var_components = 2
1016  material_field_number_of_v_var_components = ncompartments
1017  material_field_number_of_u1_var_components = 3
1018  SELECT CASE(equations_set_setup%ACTION_TYPE)
1019  CASE(equations_set_setup_start_action)
1020  equations_materials=>equations_set%MATERIALS
1021  IF(ASSOCIATED(equations_materials)) THEN
1022  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
1023  !Create the auto created materials field
1024  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
1025  & materials_field,err,error,*999)
1026  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
1027  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
1028  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1029  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
1030  & err,error,*999)
1031  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
1032  & geometric_field,err,error,*999)
1033  CALL field_number_of_variables_set(equations_materials%MATERIALS_FIELD, &
1034  & material_field_number_of_variables,err,error,*999)
1035  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type, &
1036  & field_v_variable_type,field_u1_variable_type],err,error,*999)
1037  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1038  & field_vector_dimension_type,err,error,*999)
1039  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1040  & field_dp_type,err,error,*999)
1041  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1042  & field_vector_dimension_type,err,error,*999)
1043  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1044  & field_dp_type,err,error,*999)
1045  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u1_variable_type, &
1046  & field_vector_dimension_type,err,error,*999)
1047  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u1_variable_type, &
1048  & field_dp_type,err,error,*999)
1049  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1050  & material_field_number_of_u_var_components,err,error,*999)
1051  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1052  & material_field_number_of_v_var_components,err,error,*999)
1053  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u1_variable_type, &
1054  & material_field_number_of_u1_var_components,err,error,*999)
1055  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1056  & 1,geometric_component_number,err,error,*999)
1057 
1058  !Auto-created / default is node_based_interpolation: that's an expensive default ...
1059  !Maybe default should be constant; node_based should be requested by the user \todo
1060  DO i = 1, material_field_number_of_u_var_components
1061  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1062  & i,field_node_based_interpolation,err,error,*999)
1063  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1064  & i,geometric_component_number,err,error,*999)
1065  END DO
1066  DO i = 1, material_field_number_of_v_var_components
1067  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1068  & i,field_node_based_interpolation,err,error,*999)
1069  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1070  & i,geometric_component_number,err,error,*999)
1071  END DO
1072  DO i = 1, material_field_number_of_u1_var_components
1073  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u1_variable_type, &
1074  & i,field_node_based_interpolation,err,error,*999)
1075  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u1_variable_type, &
1076  & i,geometric_component_number,err,error,*999)
1077  END DO
1078 
1079  !Default the field scaling to that of the geometric field
1080  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1081  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
1082  ELSE
1083  !Check the user specified field
1084  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
1085  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1086  CALL field_number_of_variables_check(equations_set_setup%FIELD,material_field_number_of_variables,err,error,*999)
1087  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type, &
1088  & field_v_variable_type,field_u1_variable_type],err,error,*999)
1089  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1090  & err,error,*999)
1091  CALL field_dimension_check(equations_set_setup%FIELD,field_v_variable_type,field_vector_dimension_type, &
1092  & err,error,*999)
1093  CALL field_dimension_check(equations_set_setup%FIELD,field_u1_variable_type,field_vector_dimension_type, &
1094  & err,error,*999)
1095  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1096  CALL field_data_type_check(equations_set_setup%FIELD,field_v_variable_type,field_dp_type,err,error,*999)
1097  CALL field_data_type_check(equations_set_setup%FIELD,field_u1_variable_type,field_dp_type,err,error,*999)
1098  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
1099  & material_field_number_of_u_var_components,err,error,*999)
1100  CALL field_number_of_components_check(equations_set_setup%FIELD,field_v_variable_type, &
1101  & material_field_number_of_v_var_components,err,error,*999)
1102  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u1_variable_type, &
1103  & material_field_number_of_u1_var_components,err,error,*999)
1104  ENDIF
1105  ELSE
1106  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1107  ENDIF
1108  CASE(equations_set_setup_finish_action)
1109  equations_materials=>equations_set%MATERIALS
1110  IF( ASSOCIATED(equations_materials) ) THEN
1111  IF( equations_materials%MATERIALS_FIELD_AUTO_CREATED ) THEN
1112  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
1113  !Set the default values for the materials field
1114  DO i=1,material_field_number_of_u_var_components
1115  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1116  & field_values_set_type, i, 1.0_dp, err, error, *999)
1117  ENDDO
1118  DO i=1,material_field_number_of_v_var_components
1119  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1120  & field_values_set_type, i, 0.0_dp, err, error, *999)
1121  ENDDO
1122  DO i=1,material_field_number_of_u1_var_components
1123  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u1_variable_type, &
1124  & field_values_set_type, i, 0.0_dp, err, error, *999)
1125  ENDDO
1126  ENDIF
1127  ELSE
1128  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1129  ENDIF
1130  CASE DEFAULT
1131  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1132  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1133  & " is invalid for a standard, quasistatic or ALE Darcy equation."
1134  CALL flagerror(local_error,err,error,*999)
1135  END SELECT
1136  END SELECT
1137 
1138  !-----------------------------------------------------------------
1139  ! a n a l y t i c f i e l d
1140  !-----------------------------------------------------------------
1141 
1142  CASE(equations_set_setup_analytic_type)
1143 
1144  SELECT CASE(equations_set%SPECIFICATION(3))
1145  CASE(equations_set_standard_darcy_subtype)
1146  SELECT CASE(equations_set_setup%ACTION_TYPE)
1147  !Set start action
1148  CASE(equations_set_setup_start_action)
1149  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
1150  IF(ASSOCIATED(equations_set%DEPENDENT%DEPENDENT_FIELD)) THEN
1151  IF(ASSOCIATED(equations_set%GEOMETRY%GEOMETRIC_FIELD)) THEN
1152  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1153  & number_of_dimensions,err,error,*999)
1154  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
1155  CASE(equations_set_darcy_equation_two_dim_1)
1156  !Set analytic function type
1157  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_darcy_equation_two_dim_1
1158  CASE(equations_set_darcy_equation_two_dim_2)
1159  !Set analytic function type
1160  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_darcy_equation_two_dim_2
1161  CASE(equations_set_darcy_equation_two_dim_3)
1162  !Set analytic function type
1163  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_darcy_equation_two_dim_3
1164  CASE(equations_set_darcy_equation_three_dim_1)
1165  !Set analytic function type
1166  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_darcy_equation_three_dim_1
1167  CASE(equations_set_darcy_equation_three_dim_2)
1168  !Set analytic function type
1169  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_darcy_equation_three_dim_2
1170  CASE(equations_set_darcy_equation_three_dim_3)
1171  !Set analytic function type
1172  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_darcy_equation_three_dim_3
1173  CASE(equations_set_incomp_elast_darcy_analytic_darcy)
1174  !Set analytic function type
1175  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_incomp_elast_darcy_analytic_darcy
1176  CASE DEFAULT
1177  local_error="The specified analytic function type of "// &
1178  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1179  & " is invalid for an analytic Darcy problem."
1180  CALL flagerror(local_error,err,error,*999)
1181  END SELECT
1182  ELSE
1183  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
1184  ENDIF
1185  ELSE
1186  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
1187  ENDIF
1188  ELSE
1189  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
1190  ENDIF
1191  CASE(equations_set_setup_finish_action)
1192  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
1193  IF(ASSOCIATED(equations_set%ANALYTIC%ANALYTIC_FIELD)) THEN
1194  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
1195  !--- Why finish the dependent field and not the analytic one ???
1196  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
1197  ENDIF
1198  ENDIF
1199  ELSE
1200  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
1201  ENDIF
1202  CASE DEFAULT
1203  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1204  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1205  & " is invalid for an analytic Darcy problem."
1206  CALL flagerror(local_error,err,error,*999)
1207  END SELECT
1208  CASE(equations_set_incompressible_elasticity_driven_darcy_subtype)
1209  SELECT CASE(equations_set_setup%ACTION_TYPE)
1210  !Set start action
1211  CASE(equations_set_setup_start_action)
1212  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
1213  IF(ASSOCIATED(equations_set%DEPENDENT%DEPENDENT_FIELD)) THEN
1214  IF(ASSOCIATED(equations_set%GEOMETRY%GEOMETRIC_FIELD)) THEN
1215  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1216  & number_of_dimensions,err,error,*999)
1217  !Initialise analytic parameter which stores value of time to zero - need to update this somewhere in a pre_solve routine
1218  equations_set%ANALYTIC%ANALYTIC_USER_PARAMS(1)=0.0_dp
1219  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
1220  CASE(equations_set_incomp_elast_darcy_analytic_darcy)
1221  !Set analytic function type
1222  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_incomp_elast_darcy_analytic_darcy
1223  CASE DEFAULT
1224  local_error="The specified analytic function type of "// &
1225  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1226  & " is invalid for an analytic Darcy problem."
1227  CALL flagerror(local_error,err,error,*999)
1228  END SELECT
1229  ELSE
1230  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
1231  ENDIF
1232  ELSE
1233  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
1234  ENDIF
1235  ELSE
1236  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
1237  ENDIF
1238  CASE(equations_set_setup_finish_action)
1239  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
1240  IF(ASSOCIATED(equations_set%ANALYTIC%ANALYTIC_FIELD)) THEN
1241  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
1242  !--- Why finish the dependent field and not the analytic one ???
1243  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
1244  ENDIF
1245  ENDIF
1246  ELSE
1247  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
1248  ENDIF
1249  CASE DEFAULT
1250  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1251  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1252  & " is invalid for an analytic Darcy problem."
1253  CALL flagerror(local_error,err,error,*999)
1254  END SELECT
1255  CASE DEFAULT
1256  local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1257  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1258  & " is invalid for a Darcy equation."
1259  CALL flagerror(local_error,err,error,*999)
1260  END SELECT
1261 
1262  !-----------------------------------------------------------------
1263  ! s o u r c e t y p e - include gravity at some point
1264  !-----------------------------------------------------------------
1265  CASE(equations_set_setup_source_type)
1266  SELECT CASE(equations_set%SPECIFICATION(3))
1267  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
1268  & equations_set_transient_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
1269  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
1270  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
1271  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
1272  SELECT CASE(equations_set_setup%ACTION_TYPE)
1273  CASE(equations_set_setup_start_action)
1274  equations_source=>equations_set%SOURCE
1275  IF(ASSOCIATED(equations_source)) THEN
1276  IF(equations_source%SOURCE_FIELD_AUTO_CREATED) THEN
1277  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_source% &
1278  & source_field,err,error,*999)
1279  CALL field_label_set(equations_source%SOURCE_FIELD,"Source Field",err,error,*999)
1280  CALL field_type_set_and_lock(equations_source%SOURCE_FIELD,field_general_type,err,error,*999)
1281  CALL field_dependent_type_set_and_lock(equations_source%SOURCE_FIELD,field_independent_type,err,error,*999)
1282  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1283  CALL field_mesh_decomposition_set_and_lock(equations_source%SOURCE_FIELD,geometric_decomposition, &
1284  & err,error,*999)
1285  CALL field_geometric_field_set_and_lock(equations_source%SOURCE_FIELD,equations_set%GEOMETRY% &
1286  & geometric_field,err,error,*999)
1287  CALL field_number_of_variables_set_and_lock(equations_source%SOURCE_FIELD,1,err,error,*999)
1288  CALL field_variable_types_set_and_lock(equations_source%SOURCE_FIELD,[field_u_variable_type], &
1289  & err,error,*999)
1290  CALL field_variable_label_set(equations_source%SOURCE_FIELD,field_u_variable_type,"Source", &
1291  & err,error,*999)
1292  CALL field_dimension_set_and_lock(equations_source%SOURCE_FIELD,field_u_variable_type, &
1293  & field_vector_dimension_type,err,error,*999)
1294  CALL field_data_type_set_and_lock(equations_source%SOURCE_FIELD,field_u_variable_type, &
1295  & field_dp_type,err,error,*999)
1296  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
1297  & number_of_dimensions, err, error, *999)
1298  number_of_source_components = number_of_dimensions + 1
1299  CALL field_number_of_components_set_and_lock(equations_source%SOURCE_FIELD,field_u_variable_type, &
1300  & number_of_source_components,err,error,*999)
1301 
1302  !Default the source components to the geometric interpolation setup with nodal interpolation
1303  IF(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
1304  & equations_set%SPECIFICATION(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
1305  !nodal / mesh based
1306  DO component_idx=1,number_of_dimensions !NUMBER_OF_SOURCE_COMPONENTS
1307  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1308  & component_idx,geometric_mesh_component,err,error,*999)
1309  CALL field_component_interpolation_set(equations_source%SOURCE_FIELD,field_u_variable_type, &
1310  & component_idx,field_node_based_interpolation,err,error,*999)
1311  CALL field_component_mesh_component_set(equations_source%SOURCE_FIELD,field_u_variable_type, &
1312  & component_idx,geometric_mesh_component,err,error,*999)
1313  ENDDO !component_idx
1314  !Set source component 'NUMBER_OF_DIMENSIONS + 1' according to GEOMETRIC_MESH_COMPONENT 'NUMBER_OF_DIMENSIONS'
1315  CALL field_component_interpolation_set(equations_source%SOURCE_FIELD,field_u_variable_type, &
1316  & number_of_dimensions + 1,field_node_based_interpolation,err,error,*999)
1317  CALL field_component_mesh_component_set(equations_source%SOURCE_FIELD,field_u_variable_type, &
1318  & number_of_dimensions + 1,geometric_mesh_component,err,error,*999)
1319  ENDIF
1320  !Default the field scaling to that of the geometric field
1321  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1322  CALL field_scaling_type_set(equations_source%SOURCE_FIELD,geometric_scaling_type,err,error,*999)
1323  ELSE
1324  !Check the user specified field
1325  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1326  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1327  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
1328  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
1329  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1330  & err,error,*999)
1331  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1332  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
1333  & number_of_dimensions, err, error, *999)
1334  number_of_source_components = number_of_dimensions + 1
1335  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
1336  & number_of_source_components,err,error,*999)
1337  ENDIF
1338  ELSE
1339  CALL flagerror("Equations set source is not associated.",err,error,*999)
1340  ENDIF
1341  CASE(equations_set_setup_finish_action)
1342  equations_source=>equations_set%SOURCE
1343  IF(ASSOCIATED(equations_source)) THEN
1344  IF(equations_source%SOURCE_FIELD_AUTO_CREATED) THEN
1345  !Finish creating the source field
1346  CALL field_create_finish(equations_source%SOURCE_FIELD,err,error,*999)
1347  !Set the default values for the source field
1348  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1349  & number_of_dimensions,err,error,*999)
1350  IF(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
1351  & equations_set%SPECIFICATION(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
1352  number_of_source_components = number_of_dimensions + 1
1353  ELSE
1354  number_of_source_components=0
1355  ENDIF
1356  !Now set the source values to 0.0
1357  DO component_idx=1,number_of_source_components
1358  CALL field_component_values_initialise(equations_source%SOURCE_FIELD,field_u_variable_type, &
1359  & field_values_set_type,component_idx,0.0_dp,err,error,*999)
1360  ENDDO !component_idx
1361  ENDIF
1362  ELSE
1363  CALL flagerror("Equations set source is not associated.",err,error,*999)
1364  ENDIF
1365  CASE DEFAULT
1366  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1367  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1368  & " is invalid for a standard, quasistatic or ALE Darcy equation."
1369  CALL flagerror(local_error,err,error,*999)
1370  END SELECT
1371  END SELECT
1372 
1373  !-----------------------------------------------------------------
1374  ! e q u a t i o n s t y p e
1375  !-----------------------------------------------------------------
1376  CASE(equations_set_setup_equations_type)
1377  SELECT CASE(equations_set%SPECIFICATION(3))
1378  !-----------------------------------------------------------------
1379  ! s t a t i c
1380  !-----------------------------------------------------------------
1381  CASE(equations_set_standard_darcy_subtype,equations_set_multi_compartment_darcy_subtype)
1382  SELECT CASE(equations_set_setup%ACTION_TYPE)
1383  CASE(equations_set_setup_start_action)
1384  equations_materials=>equations_set%MATERIALS
1385  IF( ASSOCIATED(equations_materials) ) THEN
1386  IF( equations_materials%MATERIALS_FINISHED ) THEN
1387  CALL equations_create_start(equations_set,equations,err,error,*999)
1388  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
1389  CALL equations_time_dependence_type_set(equations,equations_static,err,error,*999)
1390  ELSE
1391  CALL flagerror("Equations set materials has not been finished.",err,error,*999)
1392  ENDIF
1393  ELSE
1394  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1395  ENDIF
1396  CASE(equations_set_setup_finish_action)
1397  SELECT CASE(equations_set%SOLUTION_METHOD)
1398  CASE(equations_set_fem_solution_method)
1399  SELECT CASE(equations_set%SPECIFICATION(3))
1400  CASE(equations_set_multi_compartment_darcy_subtype)
1401  !!!!!THE FOLLOWING IF STATEMENT IS ILLUSTRATIVE ONLY - need to implement the equation set field thing, and make a generalised case statement
1402  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
1403  CALL equations_create_finish(equations,err,error,*999)
1404  !Create the equations mapping.
1405  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
1406  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
1407  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1408  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1409  & field_values_set_type,equations_set_field_data,err,error,*999)
1410  imy_matrix = equations_set_field_data(1)
1411  ncompartments = equations_set_field_data(2)
1412  ALLOCATE(variable_types(2*ncompartments))
1413  DO num_var=1,ncompartments
1414  variable_types(2*num_var-1)=field_u_variable_type+(field_number_of_variable_subtypes*(num_var-1))
1415  variable_types(2*num_var)=field_deludeln_variable_type+(field_number_of_variable_subtypes*(num_var-1))
1416  ENDDO
1417  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,[variable_types(2*imy_matrix-1)], &
1418  & err,error,*999)
1419  CALL equations_mapping_rhs_variable_type_set(equations_mapping,variable_types(2*imy_matrix),err,error,*999)
1420  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
1421  !Create the equations matrices
1422  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
1423  SELECT CASE(equations%SPARSITY_TYPE)
1424  CASE(equations_matrices_full_matrices)
1425  CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_block_storage_type], &
1426  & err,error,*999)
1427  CASE(equations_matrices_sparse_matrices)
1428  CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_compressed_row_storage_type], &
1429  & err,error,*999)
1430  CALL equationsmatrices_linearstructuretypeset(equations_matrices,[equations_matrix_fem_structure], &
1431  & err,error,*999)
1432  CASE DEFAULT
1433  local_error="The equations matrices sparsity type of "// &
1434  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
1435  CALL flagerror(local_error,err,error,*999)
1436  END SELECT
1437  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
1438  CASE DEFAULT
1439  !Finish the equations creation
1440  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
1441  CALL equations_create_finish(equations,err,error,*999)
1442  !Create the equations mapping.
1443  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
1444  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
1445  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,[field_u_variable_type], &
1446  & err,error,*999)
1447  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
1448  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
1449  !Create the equations matrices
1450  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
1451  SELECT CASE(equations%SPARSITY_TYPE)
1452  CASE(equations_matrices_full_matrices)
1453  CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_block_storage_type], &
1454  & err,error,*999)
1455  CASE(equations_matrices_sparse_matrices)
1456  CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_compressed_row_storage_type], &
1457  & err,error,*999)
1458  CALL equationsmatrices_linearstructuretypeset(equations_matrices,[equations_matrix_fem_structure], &
1459  & err,error,*999)
1460  CASE DEFAULT
1461  local_error="The equations matrices sparsity type of "// &
1462  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
1463  CALL flagerror(local_error,err,error,*999)
1464  END SELECT
1465  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
1466  END SELECT
1467  CASE(equations_set_bem_solution_method)
1468  CALL flagerror("Not implemented.",err,error,*999)
1469  CASE(equations_set_fd_solution_method)
1470  CALL flagerror("Not implemented.",err,error,*999)
1471  CASE(equations_set_fv_solution_method)
1472  CALL flagerror("Not implemented.",err,error,*999)
1473  CASE(equations_set_gfem_solution_method)
1474  CALL flagerror("Not implemented.",err,error,*999)
1475  CASE(equations_set_gfv_solution_method)
1476  CALL flagerror("Not implemented.",err,error,*999)
1477  CASE DEFAULT
1478  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1479  & " is invalid."
1480  CALL flagerror(local_error,err,error,*999)
1481  END SELECT
1482 
1483  CASE DEFAULT
1484  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1485  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1486  & " is invalid for a standard Darcy equation."
1487  CALL flagerror(local_error,err,error,*999)
1488  END SELECT
1489  !-----------------------------------------------------------------
1490  ! q u a s i s t a t i c and A L E
1491  !-----------------------------------------------------------------
1492  CASE(equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
1493  & equations_set_incompressible_finite_elasticity_darcy_subtype)
1494  SELECT CASE(equations_set_setup%ACTION_TYPE)
1495  CASE(equations_set_setup_start_action)
1496  equations_materials=>equations_set%MATERIALS
1497  IF( ASSOCIATED(equations_materials) ) THEN
1498  IF( equations_materials%MATERIALS_FINISHED ) THEN
1499  CALL equations_create_start(equations_set,equations,err,error,*999)
1500  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
1501  CALL equations_time_dependence_type_set(equations,equations_quasistatic,err,error,*999)
1502  ELSE
1503  CALL flagerror("Equations set materials has not been finished.",err,error,*999)
1504  ENDIF
1505  ELSE
1506  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1507  ENDIF
1508  CASE(equations_set_setup_finish_action)
1509  SELECT CASE(equations_set%SOLUTION_METHOD)
1510  CASE(equations_set_fem_solution_method)
1511  !Finish the equations creation
1512  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
1513  CALL equations_create_finish(equations,err,error,*999)
1514  !Create the equations mapping.
1515  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
1516  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
1517  SELECT CASE(equations_set%SPECIFICATION(3))
1518  CASE(equations_set_incompressible_finite_elasticity_darcy_subtype)
1519  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,[field_v_variable_type], &
1520  & err,error,*999)
1521  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_delvdeln_variable_type,err,error,*999)
1522  CASE DEFAULT
1523  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,[field_u_variable_type], &
1524  & err,error,*999)
1525  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
1526  END SELECT
1527  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
1528  !Create the equations matrices
1529  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
1530  SELECT CASE(equations%SPARSITY_TYPE)
1531  CASE(equations_matrices_full_matrices)
1532  CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_block_storage_type], &
1533  & err,error,*999)
1534  CASE(equations_matrices_sparse_matrices)
1535  CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_compressed_row_storage_type], &
1536  & err,error,*999)
1537  CALL equationsmatrices_linearstructuretypeset(equations_matrices,[equations_matrix_fem_structure], &
1538  & err,error,*999)
1539  CASE DEFAULT
1540  local_error="The equations matrices sparsity type of "// &
1541  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
1542  CALL flagerror(local_error,err,error,*999)
1543  END SELECT
1544  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
1545  CASE(equations_set_bem_solution_method)
1546  CALL flagerror("Not implemented.",err,error,*999)
1547  CASE(equations_set_fd_solution_method)
1548  CALL flagerror("Not implemented.",err,error,*999)
1549  CASE(equations_set_fv_solution_method)
1550  CALL flagerror("Not implemented.",err,error,*999)
1551  CASE(equations_set_gfem_solution_method)
1552  CALL flagerror("Not implemented.",err,error,*999)
1553  CASE(equations_set_gfv_solution_method)
1554  CALL flagerror("Not implemented.",err,error,*999)
1555  CASE DEFAULT
1556  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1557  & " is invalid."
1558  CALL flagerror(local_error,err,error,*999)
1559  END SELECT
1560  CASE DEFAULT
1561  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1562  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1563  & " is invalid for a quasistatic Darcy equation."
1564  CALL flagerror(local_error,err,error,*999)
1565  END SELECT
1566  !-----------------------------------------------------------------
1567  ! d y n a m i c
1568  !-----------------------------------------------------------------
1569  CASE(equations_set_transient_darcy_subtype,equations_set_transient_ale_darcy_subtype, &
1570  & equations_set_elasticity_darcy_inria_model_subtype,equations_set_incompressible_elasticity_driven_darcy_subtype, &
1571  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
1572  SELECT CASE(equations_set_setup%ACTION_TYPE)
1573  CASE(equations_set_setup_start_action)
1574  equations_materials=>equations_set%MATERIALS
1575  IF(ASSOCIATED(equations_materials)) THEN
1576  IF(equations_materials%MATERIALS_FINISHED) THEN
1577  CALL equations_create_start(equations_set,equations,err,error,*999)
1578  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
1579  CALL equations_time_dependence_type_set(equations,equations_first_order_dynamic,err,error,*999)
1580  ELSE
1581  CALL flagerror("Equations set materials has not been finished.",err,error,*999)
1582  ENDIF
1583  ELSE
1584  CALL flagerror("Equations materials is not associated.",err,error,*999)
1585  ENDIF
1586  CASE(equations_set_setup_finish_action)
1587  SELECT CASE(equations_set%SOLUTION_METHOD)
1588  CASE(equations_set_fem_solution_method)
1589  !Finish the equations creation
1590  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
1591  CALL equations_create_finish(equations,err,error,*999)
1592  !Create the equations mapping.
1593  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
1594  IF(equations_set%SPECIFICATION(3)==equations_set_elasticity_darcy_inria_model_subtype .OR. &
1595  & equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) THEN
1596  CALL equationsmapping_linearmatricesnumberset(equations_mapping,0,err,error,*999)
1597  ENDIF
1598  CALL equations_mapping_dynamic_matrices_set(equations_mapping,.true.,.true.,err,error,*999)
1599  SELECT CASE(equations_set%SPECIFICATION(3))
1600  CASE(equations_set_elasticity_darcy_inria_model_subtype, &
1601  & equations_set_incompressible_elasticity_driven_darcy_subtype)
1602  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_v_variable_type,err,error,*999)
1603  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_delvdeln_variable_type, &
1604  & err,error,*999)
1605  IF(equations_set%SPECIFICATION(3)==equations_set_incompressible_elasticity_driven_darcy_subtype) THEN
1606  CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
1607  ENDIF
1608  CASE(equations_set_incompressible_elast_multi_comp_darcy_subtype)
1609  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1610  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1611  & field_values_set_type,equations_set_field_data,err,error,*999)
1612  imy_matrix = equations_set_field_data(1)
1613  ncompartments = equations_set_field_data(2)
1614  CALL equationsmapping_linearmatricesnumberset(equations_mapping,ncompartments-1,err,error,*999)
1615  ALLOCATE(variable_types(2*ncompartments+2))
1616  ALLOCATE(variable_u_types(ncompartments-1))
1617  DO num_var=1,ncompartments+1
1618  variable_types(2*num_var-1)=field_u_variable_type+(field_number_of_variable_subtypes*(num_var-1))
1619  variable_types(2*num_var)=field_deludeln_variable_type+(field_number_of_variable_subtypes*(num_var-1))
1620  ENDDO
1621  num_var_count=0
1622  DO num_var=2,ncompartments+1
1623  IF((num_var-1)/=imy_matrix)THEN
1624  num_var_count=num_var_count+1
1625  variable_u_types(num_var_count)=variable_types(2*num_var-1)
1626  ENDIF
1627  ENDDO
1628  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,variable_types(2*imy_matrix+1), &
1629  & err,error,*999)
1630  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,variable_u_types,err,error,*999)
1631  CALL equations_mapping_rhs_variable_type_set(equations_mapping,variable_types(2*imy_matrix+2),err,error,*999)
1632  CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
1633  CASE DEFAULT
1634  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
1635  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type, &
1636  & err,error,*999)
1637  END SELECT
1638  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
1639  !Create the equations matrices
1640  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
1641  !Set up matrix storage and structure
1642  IF(equations%LUMPING_TYPE==equations_lumped_matrices) THEN
1643  !Set up lumping
1644  CALL equations_matrices_dynamic_lumping_type_set(equations_matrices, &
1645  & [equations_matrix_unlumped,equations_matrix_lumped],err,error,*999)
1646  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
1647  & [distributed_matrix_compressed_row_storage_type,distributed_matrix_diagonal_storage_type] &
1648  & ,err,error,*999)
1649  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
1650  & [equations_matrix_fem_structure,equations_matrix_diagonal_structure],err,error,*999)
1651  ELSE
1652  SELECT CASE(equations%SPARSITY_TYPE)
1653  CASE(equations_matrices_full_matrices)
1654  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
1655  & [distributed_matrix_block_storage_type,distributed_matrix_block_storage_type],err,error,*999)
1656  CASE(equations_matrices_sparse_matrices)
1657  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
1658  & [distributed_matrix_compressed_row_storage_type, &
1659  & distributed_matrix_compressed_row_storage_type],err,error,*999)
1660  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
1661  & [equations_matrix_fem_structure,equations_matrix_fem_structure],err,error,*999)
1662  IF(equations_set%SPECIFICATION(3)==equations_set_incompressible_elast_multi_comp_darcy_subtype)THEN
1663  ALLOCATE(coupling_matrix_storage_type(ncompartments-1))
1664  ALLOCATE(coupling_matrix_structure_type(ncompartments-1))
1665  DO num_var=1,ncompartments-1
1666  coupling_matrix_storage_type(num_var)=distributed_matrix_compressed_row_storage_type
1667  coupling_matrix_structure_type(num_var)=equations_matrix_fem_structure
1668  ENDDO
1669  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
1670  & coupling_matrix_storage_type, &
1671  & err,error,*999)
1672  CALL equationsmatrices_linearstructuretypeset(equations_matrices, &
1673  coupling_matrix_structure_type,err,error,*999)
1674  ENDIF
1675  CASE DEFAULT
1676  local_error="The equations matrices sparsity type of "// &
1677  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
1678  CALL flagerror(local_error,err,error,*999)
1679  END SELECT
1680  ENDIF
1681  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
1682  CASE DEFAULT
1683  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*", &
1684  & err,error))//" is invalid."
1685  CALL flagerror(local_error,err,error,*999)
1686  END SELECT
1687  CASE DEFAULT
1688  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1689  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1690  & " is invalid for a Darcy equation."
1691  CALL flagerror(local_error,err,error,*999)
1692  END SELECT
1693  !-----------------------------------------------------------------
1694  ! D e f a u l t
1695  !-----------------------------------------------------------------
1696  CASE DEFAULT
1697  local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1698  & " for a setup of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1699  & " is invalid for a Darcy equation."
1700  CALL flagerror(local_error,err,error,*999)
1701  END SELECT
1702 
1703  !-----------------------------------------------------------------
1704  ! c a s e d e f a u l t
1705  !-----------------------------------------------------------------
1706  CASE DEFAULT
1707  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1708  & " is invalid for a standard, quasistatic, ALE or dynamic Darcy equation."
1709  CALL flagerror(local_error,err,error,*999)
1710 
1711  END SELECT
1712  CASE DEFAULT
1713  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1714  & " does not equal a standard, quasistatic, ALE or dynamic Darcy equation subtype."
1715  CALL flagerror(local_error,err,error,*999)
1716  END SELECT
1717  ELSE
1718  CALL flagerror("Equations set is not associated.",err,error,*999)
1719  ENDIF
1720 
1721  exits("DARCY_EQUATION_EQUATIONS_SET_SETUP")
1722  RETURN
1723 999 errorsexits("DARCY_EQUATION_EQUATIONS_SET_SETUP",err,error)
1724  RETURN 1
1725  END SUBROUTINE darcy_equation_equations_set_setup
1726 
1727  !
1728  !================================================================================================================================
1729  !
1730 
1732  SUBROUTINE darcy_equation_finite_element_calculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
1733 
1734  !Argument variables
1735  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1736  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
1737  INTEGER(INTG), INTENT(OUT) :: ERR
1738  TYPE(varying_string), INTENT(OUT) :: ERROR
1739 
1740  !Local Variables
1741  INTEGER(INTG) :: FIELD_VAR_TYPE,ng,mh,mhs,mi,ms,nh,nhs,ni,ns,idxdim,num_var_count,idx_tensor
1742  INTEGER(INTG) :: my_compartment,Ncompartments,imatrix
1743  INTEGER(INTG) :: component_idx,xi_idx,derivative_idx
1744  INTEGER(INTG) :: MESH_COMPONENT_NUMBER, global_element_idx
1745  INTEGER(INTG) :: MESH_COMPONENT_1, MESH_COMPONENT_2
1746  INTEGER(INTG) :: NDOFS, NUMBER_OF_VEL_PRESS_COMPONENTS
1747  INTEGER(INTG) :: FIELD_VAR_TYPES(99)
1748  INTEGER(INTG) :: EQUATIONS_SET_SUBTYPE
1749  INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
1750  REAL(DP) :: RWG,SUM,PGMSI(3),PGNSI(3),PGM,PGN,COUPLING_PARAM
1751  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
1752  TYPE(basis_type), POINTER :: DEPENDENT_BASIS_1, DEPENDENT_BASIS_2
1753  TYPE(equations_type), POINTER :: EQUATIONS
1754  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
1755  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
1756  TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
1757  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
1758  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
1759  TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
1760  TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
1761  TYPE(equations_matrix_type), POINTER :: STIFFNESS_MATRIX, DAMPING_MATRIX
1762  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD,EQUATIONS_SET_FIELD
1763  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
1764  TYPE(decomposition_type), POINTER :: DECOMPOSITION
1765  TYPE(mesh_element_type), POINTER :: MESH_ELEMENT
1766  TYPE(boundary_conditions_variable_type), POINTER :: BOUNDARY_CONDITIONS_VARIABLE
1767  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
1768  TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
1769  TYPE(field_type), POINTER :: SOURCE_FIELD
1770  TYPE(field_variable_ptr_type) :: FIELD_VARIABLES(99)
1771  TYPE(equations_matrix_ptr_type) :: COUPLING_MATRICES(99)
1772  REAL(DP), ALLOCATABLE :: PRESSURE_COEFF(:),PRESSURE(:),GRAD_PRESSURE(:,:)
1773 
1774  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
1775  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME_1, QUADRATURE_SCHEME_2
1776  TYPE(field_interpolated_point_type), POINTER :: GEOMETRIC_INTERPOLATED_POINT,MATERIALS_INTERPOLATED_POINT
1777  TYPE(field_interpolated_point_type), POINTER :: REFERENCE_GEOMETRIC_INTERPOLATED_POINT
1778  TYPE(field_interpolation_parameters_type), POINTER :: ELASTICITY_DEPENDENT_INTERPOLATION_PARAMETERS
1779  TYPE(field_interpolated_point_type), POINTER :: ELASTICITY_DEPENDENT_INTERPOLATED_POINT
1780  TYPE(varying_string) :: LOCAL_ERROR
1781 
1782  REAL(DP):: SOURCE,INTER_COMP_SOURCE,INTER_COMP_PERM_1,INTER_COMP_PERM_2
1783  REAL(DP):: BETA_PARAM, P_SINK_PARAM
1784 
1785  REAL(DP):: PERM_OVER_VIS_PARAM, DARCY_RHO_0_F
1786  REAL(DP):: PERM_TENSOR_OVER_VIS(3,3), VIS_OVER_PERM_TENSOR(3,3), Jmat
1787  REAL(DP):: X(3), ARG(3), L, FACT
1788  REAL(DP):: LM_PRESSURE,GRAD_LM_PRESSURE(3)
1789 
1790  REAL(DP):: DXDY(3,3), DXDXI(3,3), DYDXI(3,3), DXIDY(3,3)
1791  REAL(DP):: Jxy, Jyxi
1792 
1793  REAL(DP):: Mfact, bfact, p0fact
1794 
1795  REAL(DP):: ffact !f(Jxy) of the INRIA model
1796  REAL(DP):: dfdJfact !dfdJfact = f'(Jxy) of the INRIA model
1797 
1798  REAL(DP):: C_PARAM
1799 
1800  LOGICAL :: STABILIZED
1801 
1802 
1803  !--- Parameter settings concerning the Finite Element implementation
1804  stabilized = .true.
1805  darcy%LENGTH = 10.0_dp
1806  l = darcy%LENGTH
1807 
1808  !--- testcase: default
1809  darcy%TESTCASE = 0
1810  darcy%ANALYTIC = .false.
1811 
1812  enters("DARCY_EQUATION_FINITE_ELEMENT_CALCULATE",err,error,*999)
1813 
1814  !Parameters settings for coupled elasticity Darcy INRIA model:
1815  CALL get_darcy_finite_elasticity_parameters(darcy_rho_0_f,mfact,bfact,p0fact,err,error,*999)
1816 
1817  NULLIFY(dependent_basis,geometric_basis)
1818  NULLIFY(dependent_basis_1, dependent_basis_2)
1819  NULLIFY(equations)
1820  NULLIFY(equations_mapping)
1821  NULLIFY(linear_mapping)
1822  NULLIFY(dynamic_mapping)
1823  NULLIFY(equations_matrices)
1824  NULLIFY(linear_matrices)
1825  NULLIFY(dynamic_matrices)
1826  NULLIFY(rhs_vector)
1827  NULLIFY(stiffness_matrix, damping_matrix)
1828  NULLIFY(dependent_field,geometric_field,materials_field)
1829  NULLIFY(field_variable)
1830  NULLIFY(quadrature_scheme)
1831  NULLIFY(quadrature_scheme_1, quadrature_scheme_2)
1832  NULLIFY(geometric_interpolated_point,materials_interpolated_point)
1833  NULLIFY(reference_geometric_interpolated_point)
1834  NULLIFY(elasticity_dependent_interpolation_parameters)
1835  NULLIFY(elasticity_dependent_interpolated_point)
1836  NULLIFY(decomposition,mesh_element)
1837  NULLIFY(boundary_conditions,boundary_conditions_variable)
1838  NULLIFY(source_vector,source_field)
1839  NULLIFY(equations_set_field_data)
1840 
1841  IF(ASSOCIATED(equations_set)) THEN
1842  equations=>equations_set%EQUATIONS
1843  IF(ASSOCIATED(equations)) THEN
1844  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
1845  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
1846  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
1847  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
1848  & err,error,*999)
1849  END IF
1850  equations_set_subtype=equations_set%SPECIFICATION(3)
1851  SELECT CASE(equations_set_subtype)
1852  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
1853  & equations_set_transient_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
1854  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
1855  & equations_set_incompressible_elasticity_driven_darcy_subtype,equations_set_multi_compartment_darcy_subtype, &
1856  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
1857 !!TODO: move these and scale factor adjustment out once generalised Darcy is put in.
1858  !Store all these in equations matrices/somewhere else?????
1859  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
1860  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
1861  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
1862  IF(equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
1863  & equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
1864  source_field=>equations%INTERPOLATION%SOURCE_FIELD
1865  END IF
1866 
1867  equations_matrices=>equations%EQUATIONS_MATRICES
1868  rhs_vector=>equations_matrices%RHS_VECTOR
1869  equations_mapping=>equations%EQUATIONS_MAPPING
1870 
1871  IF(equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
1872  & equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
1873  source_vector=>equations_matrices%SOURCE_VECTOR
1874  source_vector%ELEMENT_VECTOR%VECTOR = 0.0_dp
1875  END IF
1876 
1877  SELECT CASE(equations_set_subtype)
1878  CASE(equations_set_standard_darcy_subtype, equations_set_quasistatic_darcy_subtype, equations_set_ale_darcy_subtype, &
1879  equations_set_incompressible_finite_elasticity_darcy_subtype)
1880  linear_matrices=>equations_matrices%LINEAR_MATRICES
1881  stiffness_matrix=>linear_matrices%MATRICES(1)%PTR
1882 
1883  linear_mapping=>equations_mapping%LINEAR_MAPPING
1884  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1885  field_var_type=field_variable%VARIABLE_TYPE
1886 
1887  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1888 
1889  CASE(equations_set_transient_darcy_subtype,equations_set_transient_ale_darcy_subtype, &
1890  & equations_set_elasticity_darcy_inria_model_subtype, &
1891  & equations_set_incompressible_elasticity_driven_darcy_subtype)
1892  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
1893  stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
1894  damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
1895 
1896  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
1897  field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1898  field_var_type=field_variable%VARIABLE_TYPE
1899 
1900  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1901  damping_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1902 
1903  !Stuff used to check if this element is on the mesh boundary
1904  decomposition => dependent_field%DECOMPOSITION
1905  mesh_component_number = decomposition%MESH_COMPONENT_NUMBER
1906  global_element_idx = decomposition%DOMAIN(mesh_component_number)%PTR%MAPPINGS%ELEMENTS% &
1907  & local_to_global_map(element_number)
1908  mesh_element => decomposition%MESH%TOPOLOGY(mesh_component_number)%PTR%ELEMENTS%ELEMENTS(global_element_idx)
1909 
1910  CASE(equations_set_multi_compartment_darcy_subtype)
1911  equations_set_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1912  CALL field_parameter_set_data_get(equations_set_field,field_u_variable_type, &
1913  & field_values_set_type,equations_set_field_data,err,error,*999)
1914 
1915  my_compartment = equations_set_field_data(1)
1916  ncompartments = equations_set_field_data(2)
1917  linear_matrices=>equations_matrices%LINEAR_MATRICES
1918  stiffness_matrix=>linear_matrices%MATRICES(1)%PTR
1919 
1920  linear_mapping=>equations_mapping%LINEAR_MAPPING
1921  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1922  field_var_type=field_variable%VARIABLE_TYPE
1923 
1924  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1925 
1926  CASE(equations_set_elasticity_multi_compartment_darcy_inria_subtype)
1927 
1928  equations_set_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1929  CALL field_parameter_set_data_get(equations_set_field,field_u_variable_type, &
1930  & field_values_set_type,equations_set_field_data,err,error,*999)
1931 
1932  my_compartment = equations_set_field_data(1)
1933  ncompartments = equations_set_field_data(2)
1934 
1935  !if Ncompartments>99 then flag error
1936 
1937  linear_matrices=>equations_matrices%LINEAR_MATRICES
1938  linear_mapping=>equations_mapping%LINEAR_MAPPING
1939 
1940 ! DO imatrix = 1,Ncompartments
1941 ! COUPLING_MATRICES(imatrix)%PTR=>LINEAR_MATRICES%MATRICES(imatrix)%PTR
1942 ! FIELD_VARIABLES(imatrix)%PTR=>LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(imatrix)%VARIABLE
1943 ! FIELD_VAR_TYPES(imatrix)=FIELD_VARIABLES(imatrix)%PTR%VARIABLE_TYPE
1944 ! COUPLING_MATRICES(imatrix)%PTR%ELEMENT_MATRIX%MATRIX=0.0_DP
1945 ! END DO
1946 
1947  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
1948  stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
1949  damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
1950 
1951 
1952 
1953  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
1954  field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1955  field_var_type=field_variable%VARIABLE_TYPE
1956 
1957  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1958  damping_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1959  CASE(equations_set_incompressible_elast_multi_comp_darcy_subtype)
1960  equations_set_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1961  CALL field_parameter_set_data_get(equations_set_field,field_u_variable_type, &
1962  & field_values_set_type,equations_set_field_data,err,error,*999)
1963 
1964  my_compartment = equations_set_field_data(1)
1965  ncompartments = equations_set_field_data(2)
1966  !These linear matrices are actually only required if we are coupling the momentum terms too
1967  !If it is just a mass coupling, then all of the additional terms are placed in the RHS of the mass-increase equation
1968  linear_matrices=>equations_matrices%LINEAR_MATRICES
1969  linear_mapping=>equations_mapping%LINEAR_MAPPING
1970 
1971  num_var_count=0
1972  DO imatrix = 1,ncompartments
1973  IF(imatrix/=my_compartment)THEN
1974  num_var_count=num_var_count+1
1975  coupling_matrices(num_var_count)%PTR=>linear_matrices%MATRICES(num_var_count)%PTR
1976  field_variables(num_var_count)%PTR=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(num_var_count)%VARIABLE
1977  field_var_types(num_var_count)=field_variables(num_var_count)%PTR%VARIABLE_TYPE
1978  coupling_matrices(num_var_count)%PTR%ELEMENT_MATRIX%MATRIX=0.0_dp
1979  ENDIF
1980  END DO
1981 
1982  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
1983  stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
1984  damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
1985 
1986  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
1987  field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
1988  field_var_type=field_variable%VARIABLE_TYPE
1989 
1990  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1991  damping_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
1992 
1993  ALLOCATE(pressure_coeff(ncompartments))
1994 
1995  ALLOCATE(pressure(ncompartments))
1996  ALLOCATE(grad_pressure(3,ncompartments))
1997  pressure = 0.0_dp
1998  grad_pressure = 0.0_dp
1999  pressure_coeff(1)=0.25_dp
2000  pressure_coeff(2)=0.25_dp
2001  pressure_coeff(3)=0.25_dp
2002  pressure_coeff(4)=0.25_dp
2003  END SELECT
2004 
2005  !\ToDo: DEPENDENT_BASIS, DEPENDENT_BASIS_1, DEPENDENT_BASIS_2 - consistency !!!
2006 
2007  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
2008  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2009  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
2010  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2011 
2012  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
2013 
2014  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
2015  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
2016  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
2017  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
2018  IF(equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2019  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
2020  & materials_interp_parameters(field_v_variable_type)%PTR,err,error,*999)
2021  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
2022  & materials_interp_parameters(field_u1_variable_type)%PTR,err,error,*999)
2023  ENDIF
2024 
2025  IF(equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
2026  & equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2027  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
2028  & source_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
2029  END IF
2030 
2031  IF(equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
2032  & equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2033  elasticity_dependent_interpolation_parameters=>equations%INTERPOLATION% &
2034  & dependent_interp_parameters(field_u_variable_type)%PTR
2035  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
2036  & elasticity_dependent_interpolation_parameters,err,error,*999)
2037  elasticity_dependent_interpolated_point=>equations%INTERPOLATION% &
2038  & dependent_interp_point(field_u_variable_type)%PTR
2039  ENDIF
2040 
2041 
2042  SELECT CASE(equations_set_subtype)
2043  CASE(equations_set_elasticity_darcy_inria_model_subtype,equations_set_incompressible_elasticity_driven_darcy_subtype, &
2044  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
2045  number_of_vel_press_components = field_variable%NUMBER_OF_COMPONENTS - 1 !last component: mass increase
2046  CASE DEFAULT
2047  number_of_vel_press_components = field_variable%NUMBER_OF_COMPONENTS
2048  END SELECT
2049 
2050  !---------------------------------------------------------------------------------------------------------
2051  !Invoke penalty term to enforce impermeable BC
2052  ! should only be executed if THIS element lies on the surface
2053  ! (within the routine we check whether the element nodes have actually been set impermeable)
2054  SELECT CASE(equations_set_subtype)
2055  CASE(equations_set_transient_darcy_subtype,equations_set_transient_ale_darcy_subtype, &
2056  & equations_set_elasticity_darcy_inria_model_subtype, &
2057  & equations_set_incompressible_elasticity_driven_darcy_subtype)
2058  IF( mesh_element%BOUNDARY_ELEMENT ) THEN
2059  CALL darcy_equation_impermeable_bc_via_penalty(equations_set,element_number,err,error,*999)
2060  ENDIF
2061  END SELECT
2062  !---------------------------------------------------------------------------------------------------------
2063 
2064  !--- Loop over gauss points
2065  ! Given that also materials field is interpolated, ensure sufficient number of Gauss points !!!
2066  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
2067 
2068 
2069  IF(equations_set_subtype==equations_set_elasticity_darcy_inria_model_subtype.OR. &
2070  & equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
2071  & equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2072  !------------------------------------------------------------------------------
2073  !--- begin: Compute the Jacobian of the mapping
2074 
2075  !--- Interpolation of Reference Geometry
2076  CALL field_interpolation_parameters_element_get(field_initial_values_set_type,element_number, &
2077  & equations%INTERPOLATION%GEOMETRIC_INTERP_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
2078  reference_geometric_interpolated_point => equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
2079  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng, &
2080  & reference_geometric_interpolated_point,err,error,*999)
2081  !--- Retrieve local map DYDXI
2082  DO component_idx=1,dependent_basis%NUMBER_OF_XI
2083  DO xi_idx=1,dependent_basis%NUMBER_OF_XI
2084  derivative_idx=partial_derivative_first_derivative_map(xi_idx) !2,4,7
2085  dydxi(component_idx,xi_idx)=reference_geometric_interpolated_point%VALUES(component_idx,derivative_idx) !dy/dxi (y = referential)
2086  ENDDO
2087  ENDDO
2088 
2089  !--- Interpolation of (actual) Geometry and Metrics
2090  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
2091  & equations%INTERPOLATION%GEOMETRIC_INTERP_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
2092  geometric_interpolated_point => equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
2093  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng, &
2094  & geometric_interpolated_point,err,error,*999)
2095  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI, &
2096  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR,err,error,*999)
2097  !--- Retrieve local map DXDXI
2098  DO component_idx=1,dependent_basis%NUMBER_OF_XI
2099  DO xi_idx=1,dependent_basis%NUMBER_OF_XI
2100  derivative_idx=partial_derivative_first_derivative_map(xi_idx) !2,4,7
2101  dxdxi(component_idx,xi_idx)=geometric_interpolated_point%VALUES(component_idx,derivative_idx) !dx/dxi
2102  ENDDO
2103  ENDDO
2104 
2105  !--- Compute deformation gradient tensor DXDY and its Jacobian Jxy
2106  CALL invert(dydxi,dxidy,jyxi,err,error,*999) !dy/dxi -> dxi/dy
2107  CALL matrix_product(dxdxi,dxidy,dxdy,err,error,*999) !dx/dxi * dxi/dy = dx/dy (deformation gradient tensor, F)
2108  jxy=determinant(dxdy,err,error)
2109 
2110  IF( abs(jxy) < 1.0e-10_dp ) THEN
2111  local_error="DARCY_EQUATION_FINITE_ELEMENT_CALCULATE: Jacobian Jxy is smaller than 1.0E-10_DP."
2112  CALL flagerror(local_error,err,error,*999)
2113  END IF
2114 
2115  !ffact = f(Jxy) of the INRIA model, dfdJfact is not relevant here
2116  CALL evaluate_chapelle_function(jxy,ffact,dfdjfact,err,error,*999)
2117 
2118  !--- end: Compute the Jacobian of the mapping
2119  !------------------------------------------------------------------------------
2120  END IF
2121 
2122  !--- Interpolate geometric and mesh velocity field (if applicable)
2123  geometric_interpolated_point => equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
2124  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng, &
2125  & geometric_interpolated_point,err,error,*999)
2126  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI, &
2127  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR,err,error,*999)
2128 
2129  !--- Calculate 'GEOMETRIC_INTERP_PARAMETERS' from 'FIELD_VALUES_SET_TYPE'
2130  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
2131  & equations%INTERPOLATION%GEOMETRIC_INTERP_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
2132  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng, &
2133  & geometric_interpolated_point,err,error,*999)
2134 
2135 
2136 ! !--- Material Settings ---!
2137 ! !*** If material is variable, need to account for this in deriving the variational statement ***!
2138 
2139 
2140  !--- Interpolate materials field
2141  !Get the Darcy permeability
2142  materials_interpolated_point => equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR
2143  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng, &
2144  & materials_interpolated_point,err,error,*999)
2145  IF(equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2146  !Get the intercompartmental permeabilities
2147  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
2148  & materials_interp_point(field_v_variable_type)%PTR,err,error,*999)
2149  !Get the material parameters for the constitutive law for each Darcy compartment (for determining the partial pressures)
2150  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
2151  & materials_interp_point(field_u1_variable_type)%PTR,err,error,*999)
2152  ENDIF
2153 
2154  IF(equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
2155  & equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2156  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
2157  & source_interp_point(field_u_variable_type)%PTR,err,error,*999)
2158  END IF
2159 
2160  SELECT CASE(equations_set_subtype)
2161  CASE(equations_set_standard_darcy_subtype,equations_set_quasistatic_darcy_subtype,equations_set_ale_darcy_subtype)
2162  !scalar permeability/viscosity
2163  perm_tensor_over_vis=0.0_dp
2164  perm_tensor_over_vis(1,1) = materials_interpolated_point%VALUES(2,no_part_deriv)
2165  perm_tensor_over_vis(2,2) = materials_interpolated_point%VALUES(2,no_part_deriv)
2166  perm_tensor_over_vis(3,3) = materials_interpolated_point%VALUES(2,no_part_deriv)
2167  !Multiply by porosity
2168  perm_tensor_over_vis=perm_tensor_over_vis*materials_interpolated_point%VALUES(1,no_part_deriv)
2169  CASE DEFAULT
2170  !symmetric permeability/viscosity tensor
2171  perm_tensor_over_vis(1,1) = materials_interpolated_point%VALUES(2,no_part_deriv)
2172  perm_tensor_over_vis(1,2) = materials_interpolated_point%VALUES(3,no_part_deriv)
2173  perm_tensor_over_vis(1,3) = materials_interpolated_point%VALUES(4,no_part_deriv)
2174  perm_tensor_over_vis(2,2) = materials_interpolated_point%VALUES(5,no_part_deriv)
2175  perm_tensor_over_vis(2,3) = materials_interpolated_point%VALUES(6,no_part_deriv)
2176  perm_tensor_over_vis(3,3) = materials_interpolated_point%VALUES(7,no_part_deriv)
2177 
2178  perm_tensor_over_vis(2,1) = perm_tensor_over_vis(1,2)
2179  perm_tensor_over_vis(3,1) = perm_tensor_over_vis(1,3)
2180  perm_tensor_over_vis(3,2) = perm_tensor_over_vis(2,3)
2181  END SELECT
2182 
2183  IF(diagnostics3) THEN
2184  IF(idebug2) THEN
2185  CALL write_string_value(diagnostic_output_type,"MATERIALS_INTERPOLATED_POINT%VALUES(1,NO_PART_DERIV) = ", &
2186  & materials_interpolated_point%VALUES(1,no_part_deriv),err,error,*999)
2187  CALL write_string_value(diagnostic_output_type,"MATERIALS_INTERPOLATED_POINT%VALUES(2,NO_PART_DERIV) = ", &
2188  & materials_interpolated_point%VALUES(2,no_part_deriv),err,error,*999)
2189  CALL write_string(diagnostic_output_type,"",err,error,*999)
2190  idebug2 = .false.
2191  ENDIF
2192  ENDIF
2193 
2194  jmat = determinant(perm_tensor_over_vis,err,error)
2195  IF(jmat>zero_tolerance) THEN
2196  CALL invert(perm_tensor_over_vis,vis_over_perm_tensor,jmat,err,error,*999)
2197  ELSE
2198  vis_over_perm_tensor = 0.0_dp
2199  DO idx_tensor=1,3
2200  vis_over_perm_tensor(idx_tensor,idx_tensor) = 1.0e10_dp
2201  END DO
2202 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE, &
2203 ! & "WARNING: Jmat<ZERO_TOLERANCE - Thus setting VIS_OVER_PERM_TENSOR(i,i) = 1.0e10_DP",ERR,ERROR,*999)
2204  END IF
2205 
2206 
2207  !Two parameters that are used only for TESTCASE==3: VenousCompartment problem: Exclude this, too specific ???
2208  beta_param = - darcy%PERM_OVER_VIS * (2.0_dp * pi / darcy%LENGTH) * (2.0_dp * pi / darcy%LENGTH)
2209  p_sink_param = darcy%P_SINK
2210 
2211 
2212  IF(equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype) THEN
2213  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng, &
2214  & elasticity_dependent_interpolated_point,err,error,*999)
2215  !Mind the sign !!!
2216  !The minus sign derives from the convention of using "+ P * Jznu * AZU(i,j)"
2217  ! in the constitutive law in FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR
2218  lm_pressure = -elasticity_dependent_interpolated_point%VALUES(4,no_part_deriv)
2219  DO xi_idx=1,dependent_basis%NUMBER_OF_XI
2220  derivative_idx=partial_derivative_first_derivative_map(xi_idx) !2,4,7
2221  !gradient wrt. element coordinates xi
2222  grad_lm_pressure(xi_idx) = -elasticity_dependent_interpolated_point%VALUES(4,derivative_idx)
2223  ENDDO
2224  ENDIF
2225 
2226  !For multi-compartment model - determine pressure from partial derivative of constitutive law
2227  IF(equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2228  write(*,*) 'NEED CONSTITUTIVE LAWS HERE!!!! THE FOLLOWING IS PLACEHOLDER ONLY!'
2229  !BEGIN PLACEHOLDER
2230  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng, &
2231  & elasticity_dependent_interpolated_point,err,error,*999)
2232  !Mind the sign !!!
2233  lm_pressure = -elasticity_dependent_interpolated_point%VALUES(4,no_part_deriv)
2234  DO xi_idx=1,dependent_basis%NUMBER_OF_XI
2235  derivative_idx=partial_derivative_first_derivative_map(xi_idx) !2,4,7
2236  !gradient wrt. element coordinates xi
2237  grad_lm_pressure(xi_idx) = -elasticity_dependent_interpolated_point%VALUES(4,derivative_idx)
2238  ENDDO
2239  !loop over compartments to determine the pressure in each one - this could be quite inefficient, as it will be calculated several times over
2240  !unless calculate the pressures in a pre-solve and store them in extra components/variables of the dependent field
2241  !these pressures should really be known immediately after the finite elasticity solve and not determined here
2242  !END PLACEHOLDER
2243  !The following pressure_coeff matrix is just for testing purposes and ultimately will be replaced with functions and materials field parameters (for present, sum of
2244  !coefficients should be 1).
2245 
2246  DO imatrix=1,ncompartments
2247 
2248  pressure(imatrix) = pressure_coeff(imatrix)*lm_pressure
2249  DO xi_idx=1,dependent_basis%NUMBER_OF_XI
2250  grad_pressure(xi_idx,imatrix) = pressure_coeff(imatrix)*grad_lm_pressure(xi_idx)
2251  ENDDO
2252  ENDDO
2253  ENDIF
2254 
2255 
2256 ! RWG = EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN * QUADRATURE_SCHEME%GAUSS_WEIGHTS(ng)
2257 
2258  !Loop over element rows
2259  mhs=0
2260  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
2261 
2262  mesh_component_1 = field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
2263  dependent_basis_1 => dependent_field%DECOMPOSITION%DOMAIN(mesh_component_1)%PTR% &
2264  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2265  quadrature_scheme_1 => dependent_basis_1%QUADRATURE% &
2266  & quadrature_scheme_map(basis_default_quadrature_scheme)%PTR
2267  rwg = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN * &
2268  & quadrature_scheme_1%GAUSS_WEIGHTS(ng)
2269 
2270  DO ms=1,dependent_basis_1%NUMBER_OF_ELEMENT_PARAMETERS
2271  mhs=mhs+1
2272 
2273  !===================================================================================================================
2274  !STIFFNESS_MATRIX
2275  IF(stiffness_matrix%UPDATE_MATRIX) THEN
2276 
2277  !Loop over element columns
2278  nhs=0
2279  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
2280 
2281  mesh_component_2 = field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
2282  dependent_basis_2 => dependent_field%DECOMPOSITION%DOMAIN(mesh_component_2)%PTR% &
2283  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2284  !--- We cannot use two different quadrature schemes here !!!
2285  quadrature_scheme_2 => dependent_basis_2%QUADRATURE% &
2286  & quadrature_scheme_map(basis_default_quadrature_scheme)%PTR
2287  !RWG = EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN * &
2288  ! & QUADRATURE_SCHEME_2%GAUSS_WEIGHTS(ng)
2289 
2290  DO ns=1,dependent_basis_2%NUMBER_OF_ELEMENT_PARAMETERS
2291  nhs=nhs+1
2292 
2293  SELECT CASE(equations_set_subtype)
2294  !====================================================================================================
2295  ! i n c o m p r e s s i b l e e l a s t i c i t y d r i v e n D a r c y : M A T R I C E S
2296  CASE(equations_set_incompressible_elasticity_driven_darcy_subtype)
2297  !-------------------------------------------------------------------------------------------------------------
2298  !velocity test function, velocity trial function
2299  IF(mh==nh.AND.nh<field_variable%NUMBER_OF_COMPONENTS) THEN
2300 
2301  sum = 0.0_dp
2302 
2303  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2304  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2305 
2306  sum = sum + vis_over_perm_tensor( mh, nh ) * pgm * pgn
2307  !MIND: double check the matrix index order: (mh, nh) or (nh, mh)
2308  !within this conditional: mh==nh anyway
2309 
2310  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2311  & sum * rwg
2312  !-------------------------------------------------------------------------------------------------------------
2313  !mass-increase test function, velocity trial function
2314  ELSE IF(mh==field_variable%NUMBER_OF_COMPONENTS.AND.nh<field_variable%NUMBER_OF_COMPONENTS) THEN
2315 
2316  sum = 0.0_dp
2317 
2318  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2319  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2320 
2321  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2322  pgmsi(mi)=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(mi),ng)
2323  ENDDO !mi
2324 
2325  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2326  pgnsi(ni)=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
2327  ENDDO !ni
2328 
2329  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2330  sum = sum + pgm * pgnsi(ni) * &
2331  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nh)
2332  ENDDO !ni
2333 
2334  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2335  & sum * rwg
2336  ENDIF
2337  !=====================================================================================================================
2338  !DAMPING_MATRIX
2339  IF(damping_matrix%UPDATE_MATRIX) THEN
2340  !MASS-INCREASE test function, mass-increase trial function
2341  IF(mh==nh.AND.nh==field_variable%NUMBER_OF_COMPONENTS) THEN
2342  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2343  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2344 
2345  sum = 0.0_dp
2346 
2347  !To integrate the mass-increase term in the reference configuration, we divide by Jxy.
2348  sum = pgm * pgn / (jxy * darcy_rho_0_f)
2349 
2350  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2351  & sum * rwg
2352  END IF
2353 
2354 ! !Try out adding the inertia term ...
2355 ! IF(mh==nh.AND.mh<FIELD_VARIABLE%NUMBER_OF_COMPONENTS) THEN
2356 ! PGM=QUADRATURE_SCHEME_1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng)
2357 ! PGN=QUADRATURE_SCHEME_2%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng)
2358 !
2359 ! SUM = 0.0_DP
2360 !
2361 ! SUM = PGM*PGN*DARCY_RHO_0_F
2362 !
2363 ! DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs) = DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2364 ! & SUM * RWG
2365 ! END IF
2366 
2367  END IF
2368 
2369  ! matrices for multi-compartment poroelastic equations
2370  CASE(equations_set_incompressible_elast_multi_comp_darcy_subtype)
2371  !velocity test function, velocity trial function
2372  IF(mh==nh.AND.nh<field_variable%NUMBER_OF_COMPONENTS) THEN
2373 
2374  sum = 0.0_dp
2375 
2376  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2377  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2378 
2379  sum = sum + vis_over_perm_tensor( mh, nh ) * pgm * pgn
2380  !MIND: double check the matrix index order: (mh, nh) or (nh, mh)
2381  !within this conditional: mh==nh anyway
2382 
2383  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2384  & sum * rwg
2385  !-------------------------------------------------------------------------------------------------------------
2386  !mass-increase test function, velocity trial function
2387  ELSE IF(mh==field_variable%NUMBER_OF_COMPONENTS.AND.nh<field_variable%NUMBER_OF_COMPONENTS) THEN
2388 
2389  sum = 0.0_dp
2390 
2391  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2392  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2393 
2394  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2395  pgmsi(mi)=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(mi),ng)
2396  ENDDO !mi
2397 
2398  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2399  pgnsi(ni)=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
2400  ENDDO !ni
2401 
2402  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2403  sum = sum + pgm * pgnsi(ni) * &
2404  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nh)
2405  ENDDO !ni
2406 
2407  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2408  & sum * rwg
2409  ENDIF
2410  !=====================================================================================================================
2411  !DAMPING_MATRIX
2412  IF(damping_matrix%UPDATE_MATRIX) THEN
2413  !MASS-INCREASE test function, mass-increase trial function
2414  IF(mh==nh.AND.nh==field_variable%NUMBER_OF_COMPONENTS) THEN
2415  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2416  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2417 
2418  sum = 0.0_dp
2419 
2420  !To integrate the mass-increase term in the reference configuration, we divide by Jxy.
2421  sum = pgm * pgn / (jxy * darcy_rho_0_f)
2422 
2423  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2424  & sum * rwg
2425  END IF
2426 
2427 ! !Try out adding the inertia term ...
2428 ! IF(mh==nh.AND.mh<FIELD_VARIABLE%NUMBER_OF_COMPONENTS) THEN
2429 ! PGM=QUADRATURE_SCHEME_1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng)
2430 ! PGN=QUADRATURE_SCHEME_2%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng)
2431 !
2432 ! SUM = 0.0_DP
2433 !
2434 ! SUM = PGM*PGN*DARCY_RHO_0_F
2435 !
2436 ! DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs) = DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2437 ! & SUM * RWG
2438 ! END IF
2439 
2440  END IF
2441 
2442 
2443  !=================================================================================
2444  ! d e f a u l t : M A T R I C E S
2445  CASE DEFAULT
2446  !-------------------------------------------------------------------------------------------------------------
2447  !velocity test function, velocity trial function
2448  IF(mh==nh.AND.nh<number_of_vel_press_components) THEN
2449 
2450  sum = 0.0_dp
2451 
2452  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2453  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2454 
2455  sum = sum + vis_over_perm_tensor( mh, nh ) * pgm * pgn
2456  !MIND: double check the matrix index order: (mh, nh) or (nh, mh)
2457  !within this conditional: mh==nh anyway
2458 
2459  IF( stabilized ) THEN
2460  sum = sum - 0.5_dp * vis_over_perm_tensor( mh, nh ) * pgm * pgn
2461  END IF
2462 
2463  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2464  & sum * rwg
2465 
2466  !-------------------------------------------------------------------------------------------------------------
2467  !velocity test function, pressure trial function
2468  ELSE IF(mh<number_of_vel_press_components.AND.nh==number_of_vel_press_components) THEN
2469 
2470  sum = 0.0_dp
2471 
2472  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2473  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2474 
2475  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2476  pgmsi(mi)=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(mi),ng)
2477  ENDDO !mi
2478 
2479  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2480  pgnsi(ni)=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
2481  ENDDO !ni
2482 
2483  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2484  sum = sum - pgmsi(mi) * pgn * &
2485  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(mi,mh)
2486  ENDDO !mi
2487 
2488  IF( stabilized ) THEN
2489  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2490  sum = sum - 0.5_dp * pgm * pgnsi(ni) * &
2491  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,mh)
2492  ENDDO !ni
2493  END IF
2494 
2495  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2496  & sum * rwg
2497 
2498  !-------------------------------------------------------------------------------------------------------------
2499  !pressure test function, velocity trial function
2500  ELSE IF(mh==number_of_vel_press_components.AND.nh<number_of_vel_press_components) THEN
2501 
2502  sum = 0.0_dp
2503 
2504  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2505  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2506 
2507  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2508  pgmsi(mi)=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(mi),ng)
2509  ENDDO !mi
2510 
2511  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2512  pgnsi(ni)=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
2513  ENDDO !ni
2514 
2515  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2516  sum = sum + pgm * pgnsi(ni) * &
2517  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nh)
2518  ENDDO !ni
2519 
2520  IF( stabilized ) THEN
2521  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2522  sum = sum + 0.5_dp * pgmsi(mi) * pgn * &
2523  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(mi,nh)
2524  ENDDO !mi
2525  END IF
2526 
2527  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2528  & sum * rwg
2529 
2530  !-------------------------------------------------------------------------------------------------------------
2531  !pressure test function, pressure trial function
2532  ELSE IF(mh==nh.AND.nh==number_of_vel_press_components) THEN
2533 
2534  sum = 0.0_dp
2535 
2536  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2537  pgmsi(mi)=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(mi),ng)
2538  ENDDO !mi
2539 
2540  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2541  pgnsi(ni)=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
2542  ENDDO !ni
2543 
2544  IF( stabilized ) THEN
2545  DO idxdim =1,dependent_basis_1%NUMBER_OF_XI !number space dimension equiv. number of xi
2546  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2547  DO ni=1,dependent_basis_2%NUMBER_OF_XI
2548  sum = sum + 0.5_dp * perm_tensor_over_vis( idxdim, idxdim ) * pgmsi(mi) * pgnsi(ni) * &
2549  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR% &
2550  & dxi_dx(mi,idxdim) * equations%INTERPOLATION% &
2551  & geometric_interp_point_metrics(field_u_variable_type)%PTR%DXI_DX(ni,idxdim)
2552  ENDDO !ni
2553  ENDDO !mi
2554  ENDDO !idxdim
2555  END IF
2556 
2557  IF( darcy%TESTCASE == 3 ) THEN
2558  !This forms part of the pressure-dependent source term,
2559  !thus it enters the LHS
2560  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2561  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2562 
2563  sum = sum + beta_param * pgm * pgn
2564  END IF
2565 
2566  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2567  & sum * rwg
2568 
2569  !-------------------------------------------------------------------------------------------------------------
2570  !For the INRIA model, and: mass-increase test function, pressure trial function
2571  ELSE IF(equations_set_subtype==equations_set_elasticity_darcy_inria_model_subtype.AND. &
2572  & mh==field_variable%NUMBER_OF_COMPONENTS.AND.nh==number_of_vel_press_components) THEN
2573 
2574  sum = 0.0_dp
2575 
2576  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2577  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2578 
2579  sum = sum - pgm * pgn / (mfact * ffact)
2580 
2581  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2582  & sum * rwg
2583 
2584  !-------------------------------------------------------------------------------------------------------------
2585  !For the INRIA model, and: mass-increase test function, mass-increase trial function
2586  ELSE IF(equations_set_subtype==equations_set_elasticity_darcy_inria_model_subtype.AND. &
2587  & mh==nh.AND.nh==field_variable%NUMBER_OF_COMPONENTS) THEN
2588 
2589  sum = 0.0_dp
2590 
2591  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2592  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2593 
2594  sum = sum + pgm * pgn / darcy_rho_0_f
2595 
2596  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2597  & sum * rwg
2598 
2599  !-------------------------------------------------------------------------------------------------------------
2600  ELSE
2601 
2602  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = 0.0_dp
2603 
2604  ENDIF
2605 
2606  !=====================================================================================================================
2607  ! DAMPING_MATRIX
2608 
2609  IF(equations_set_subtype==equations_set_transient_darcy_subtype) THEN
2610  IF(damping_matrix%UPDATE_MATRIX) THEN
2611  IF(mh==nh.AND.mh<number_of_vel_press_components) THEN
2612  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2613  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2614 
2615  sum = 0.0_dp
2616 
2617  sum = pgm*pgn*darcy_rho_0_f
2618 
2619  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2620  & sum * rwg
2621  END IF
2622  END IF
2623  ELSE IF(equations_set_subtype==equations_set_elasticity_darcy_inria_model_subtype) THEN
2624  IF(damping_matrix%UPDATE_MATRIX) THEN
2625  !pressure test function, mass-increase trial function
2626  IF(mh==number_of_vel_press_components.AND.nh==field_variable%NUMBER_OF_COMPONENTS) THEN
2627  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2628  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2629 
2630  sum = 0.0_dp
2631 
2632  sum = pgm * pgn / (jxy * darcy_rho_0_f)
2633 
2634  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) = damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2635  & sum * rwg
2636  END IF
2637  END IF
2638  END IF
2639 
2640  END SELECT
2641  ! e n d s e l e c t EQUATIONS_SET_SUBTYPE
2642  !=================================================================================
2643 
2644  ENDDO !ns
2645  ENDDO !nh
2646  ENDIF
2647  !===================================================================================================================
2648  !RHS_VECTOR
2649  IF(rhs_vector%UPDATE_VECTOR) THEN
2650 
2651  SELECT CASE(equations_set_subtype)
2652  !==========================================================================================
2653  ! i n c o m p r e s s i b l e e l a s t i c i t y d r i v e n D a r c y : R H S
2654  CASE(equations_set_incompressible_elasticity_driven_darcy_subtype)
2655 
2656  !-----------------------------------------------------------------------------------------------------------------
2657  !velocity test function
2658  IF( mh<field_variable%NUMBER_OF_COMPONENTS ) THEN
2659 
2660  sum = 0.0_dp
2661 
2662  pgm = quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2663 
2664  !Term arising from the pressure / Lagrange Multiplier of elasticity (given):
2665  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2666  sum = sum - pgm * grad_lm_pressure(mi) * &
2667  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(mi,mh)
2668  ENDDO !mi
2669 
2670  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2671 
2672  !-----------------------------------------------------------------------------------------------------------------
2673  !mass-increase test function
2674  ELSE IF( mh==field_variable%NUMBER_OF_COMPONENTS ) THEN
2675 
2676  sum = 0.0_dp
2677 
2678  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2679 
2680  ! + possible SOURCE AND SINK TERMS
2681  source = 0.0_dp
2682 
2683  sum = sum + pgm * source
2684 
2685  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2686 
2687  ELSE
2688 
2689  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = 0.0_dp
2690 
2691  END IF
2692  !-------------------------------------------------------------------------------------------------------------
2693  CASE(equations_set_incompressible_elast_multi_comp_darcy_subtype)
2694  !-----------------------------------------------------------------------------------------------------------------
2695  !velocity test function
2696  IF( mh<field_variable%NUMBER_OF_COMPONENTS ) THEN
2697 
2698  sum = 0.0_dp
2699 
2700  pgm = quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2701 
2702  !Term arising from the pressure / Lagrange Multiplier of elasticity (given):
2703  !TO DO- need to read different grad p depending on the compartment of interest
2704  DO mi=1,dependent_basis_1%NUMBER_OF_XI
2705 ! SUM = SUM - PGM * GRAD_LM_PRESSURE(mi) * &
2706 ! & EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(FIELD_U_VARIABLE_TYPE)%PTR%DXI_DX(mi,mh)
2707  !this is the pressure gradient for the appropriate compartment
2708  sum = sum - pgm * grad_pressure(mi,my_compartment) * &
2709  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(mi,mh)
2710  ENDDO !mi
2711 
2712  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2713 
2714  !-----------------------------------------------------------------------------------------------------------------
2715  !mass-increase test function
2716  ELSE IF( mh==field_variable%NUMBER_OF_COMPONENTS ) THEN
2717 
2718  sum = 0.0_dp
2719 
2720  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2721 
2722  ! n o s o u r c e
2723  !source terms need to be converted to use source field & vector
2724  source = 0.0_dp
2725 
2726 
2727  !Add in the source/sink terms due to the pressure difference between compartments
2728  DO imatrix=1,ncompartments
2729  IF(imatrix/=my_compartment) THEN
2730  !Interpolate the coupling material parameter from the V variable type of the materials field
2731  inter_comp_perm_1=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_v_variable_type)%PTR% &
2732  & values(my_compartment,no_part_deriv)
2733  inter_comp_perm_2=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_v_variable_type)%PTR% &
2734  & values(imatrix,no_part_deriv)
2735  !Source term is coefficient*(p(my_compartment) - p(imatrix))
2736  inter_comp_source=-inter_comp_perm_1*pressure(my_compartment) + inter_comp_perm_2*pressure(imatrix)
2737  ENDIF
2738  ENDDO
2739 
2740 
2741  sum = sum + pgm * source + pgm * inter_comp_source
2742 
2743  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2744 
2745  ELSE
2746 
2747  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = 0.0_dp
2748 
2749  END IF
2750  !=================================================================================
2751  ! d e f a u l t : R H S
2752  CASE DEFAULT
2753  !-----------------------------------------------------------------------------------------------------------------
2754  !velocity test function
2755  IF( mh<number_of_vel_press_components ) THEN
2756 
2757  sum = 0.0_dp
2758 
2759  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2760 
2761  !-----------------------------------------------------------------------------------------------------------------
2762  !pressure test function
2763  ELSE IF( mh==number_of_vel_press_components ) THEN
2764 
2765  sum = 0.0_dp
2766 
2767  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2768 
2769  ! n o s o u r c e
2770  source = 0.0_dp
2771 
2772  sum = sum + pgm * source
2773 
2774  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2775 
2776  !-------------------------------------------------------------------------------------------------------------
2777  !For the INRIA model, and: mass-increase test function
2778  ELSE IF(equations_set_subtype==equations_set_elasticity_darcy_inria_model_subtype.AND. &
2779  & mh==field_variable%NUMBER_OF_COMPONENTS) THEN
2780 
2781  sum = 0.0_dp
2782 
2783  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2784 
2785  sum = sum - pgm * bfact * (1.0_dp - jxy)
2786 
2787  sum = sum - pgm * p0fact / (mfact * ffact)
2788 
2789  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2790 
2791  ELSE
2792 
2793  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs) = 0.0_dp
2794 
2795  END IF
2796  !-------------------------------------------------------------------------------------------------------------
2797  END SELECT
2798  ! e n d s e l e c t EQUATIONS_SET_SUBTYPE
2799  !=================================================================================
2800 
2801  END IF
2802 
2803  IF(ASSOCIATED(source_vector)) THEN
2804  IF(source_vector%UPDATE_VECTOR) THEN
2805  IF(equations_set_subtype==equations_set_incompressible_elasticity_driven_darcy_subtype .OR. &
2806  & equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2807 
2808  c_param=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(mh, no_part_deriv)
2809 
2810  !IF(ABS(C_PARAM)>1.0E-08) WRITE(*,*)'C_PARAM = ',C_PARAM
2811 
2812  sum = 0.0_dp
2813  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2814  sum = sum + pgm * c_param
2815  source_vector%ELEMENT_VECTOR%VECTOR(mhs) = source_vector%ELEMENT_VECTOR%VECTOR(mhs) + sum * rwg
2816  ENDIF
2817  END IF
2818  END IF
2819  ENDDO !ms
2820  ENDDO !mh
2821 
2822  IF(equations_set_subtype==equations_set_incompressible_elast_multi_comp_darcy_subtype) THEN
2823  !Calculate the momentum coupling matrices
2824 
2825  !Loop over element rows
2826  mhs=0
2827  DO mh=1,field_variable%NUMBER_OF_COMPONENTS !field_variable is the variable associated with the equations set under consideration
2828 
2829  mesh_component_1 = field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
2830  dependent_basis_1 => dependent_field%DECOMPOSITION%DOMAIN(mesh_component_1)%PTR% &
2831  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2832  quadrature_scheme_1 => dependent_basis_1%QUADRATURE% &
2833  & quadrature_scheme_map(basis_default_quadrature_scheme)%PTR
2834  rwg = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN * &
2835  & quadrature_scheme_1%GAUSS_WEIGHTS(ng)
2836 
2837  DO ms=1,dependent_basis_1%NUMBER_OF_ELEMENT_PARAMETERS
2838  mhs=mhs+1
2839 
2840  num_var_count=0
2841  DO imatrix = 1,ncompartments
2842  IF(imatrix/=my_compartment)THEN
2843  num_var_count=num_var_count+1
2844 
2845 !need to test for the case where imatrix==mycompartment
2846 !the coupling terms then needs to be added into the stiffness matrix
2847  IF(coupling_matrices(num_var_count)%PTR%UPDATE_MATRIX) THEN
2848 
2849  !Loop over element columns
2850  nhs=0
2851  DO nh=1,field_variables(num_var_count)%PTR%NUMBER_OF_COMPONENTS
2852 
2853  mesh_component_2 = field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
2854  dependent_basis_2 => dependent_field%DECOMPOSITION%DOMAIN(mesh_component_2)%PTR% &
2855  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2856  !--- We cannot use two different quadrature schemes here !!!
2857  quadrature_scheme_2 => dependent_basis_2%QUADRATURE% &
2858  & quadrature_scheme_map(basis_default_quadrature_scheme)%PTR
2859  !RWG = EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN * &
2860  ! & QUADRATURE_SCHEME_2%GAUSS_WEIGHTS(ng)
2861 
2862  DO ns=1,dependent_basis_2%NUMBER_OF_ELEMENT_PARAMETERS
2863  nhs=nhs+1
2864 
2865 ! !-------------------------------------------------------------------------------------------------------------
2866 ! !concentration test function, concentration trial function
2867 ! !For now, this is only a dummy implementation - this still has to be properly set up.
2868 ! IF(mh==nh.AND.nh<NUMBER_OF_VEL_PRESS_COMPONENTS) THEN ! don't need this for diffusion equation
2869 
2870 ! SUM = 0.0_DP
2871 
2872  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2873  pgn=quadrature_scheme_2%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
2874 
2875  !Get the coupling coefficients
2876  coupling_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_v_variable_type)%PTR% &
2877  & values(imatrix,no_part_deriv)
2878 
2879 ! SUM = SUM + COUPLING_PARAM * PGM * PGN
2880 
2881  coupling_matrices(num_var_count)%PTR%ELEMENT_MATRIX%MATRIX(mhs,nhs) = &
2882  & coupling_matrices(num_var_count)%PTR%ELEMENT_MATRIX%MATRIX(mhs,nhs) + &
2883  & coupling_param * pgm * pgn * rwg
2884 ! ENDIF
2885 
2886  ENDDO !ns
2887  ENDDO !nh
2888  ENDIF
2889  ENDIF
2890  ENDDO !imatrix
2891  ENDDO !ms
2892  ENDDO !mh
2893  ENDIF
2894 
2895 
2896  !-----------------------------------------------------------------------------------------------------------------------------------
2897  ! RIGHT HAND SIDE FOR ANALYTIC SOLUTION
2898  !-----------------------------------------------------------------------------------------------------------------------------------
2899 
2900  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
2901  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_two_dim_1.OR. &
2902  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_two_dim_2.OR. &
2903  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_two_dim_3.OR. &
2904  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_three_dim_1.OR. &
2905  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_three_dim_2.OR. &
2906  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_three_dim_3) THEN
2907 
2908  mhs=0
2909  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
2910  mesh_component_1=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
2911  dependent_basis_1=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_1)%PTR% &
2912  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2913  quadrature_scheme_1=>dependent_basis_1%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
2914  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
2915  & quadrature_scheme_1%GAUSS_WEIGHTS(ng)
2916  DO ms=1,dependent_basis_1%NUMBER_OF_ELEMENT_PARAMETERS
2917  mhs=mhs+1
2918  pgm=quadrature_scheme_1%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
2919  !note mh value derivative
2920  sum=0.0_dp
2921 
2922  x(1) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,1)
2923  x(2) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,1)
2924  IF(dependent_basis_1%NUMBER_OF_XI==3) THEN
2925  x(3) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,1)
2926  END IF
2927  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_two_dim_1) THEN
2928  sum=0.0_dp
2929  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_two_dim_2) THEN
2930  IF(mh==3) THEN
2931  fact = perm_over_vis_param / l
2932  arg(1) = x(1) / l
2933  arg(2) = x(2) / l
2934  source = -2.0_dp / l * fact * exp( arg(1) ) * exp( arg(2) )
2935  sum = pgm * source
2936  ELSE
2937  sum = 0.0_dp
2938  ENDIF
2939  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_two_dim_3) THEN
2940  IF(mh==3) THEN
2941  fact = 2.0_dp * pi * perm_over_vis_param / l
2942  arg(1) = 2.0_dp * pi * x(1) / l
2943  arg(2) = 2.0_dp * pi * x(2) / l
2944  source = +2.0_dp * (2.0_dp * pi / l) * fact * sin( arg(1) ) * sin( arg(2) )
2945  sum = pgm * source
2946  ELSE
2947  sum = 0.0_dp
2948  ENDIF
2949  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_three_dim_1) THEN
2950  sum=0.0_dp
2951  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_three_dim_2) THEN
2952  IF(mh==4) THEN
2953  fact = perm_over_vis_param / l
2954  arg(1) = x(1) / l
2955  arg(2) = x(2) / l
2956  arg(3) = x(3) / l
2957  source = -3.0_dp / l * fact * exp( arg(1) ) * exp( arg(2) ) * exp( arg(3) )
2958  sum = pgm * source
2959  ELSE
2960  sum = 0.0_dp
2961  ENDIF
2962  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_darcy_equation_three_dim_3) THEN
2963  IF(mh==4) THEN
2964  fact = 2.0_dp * pi * perm_over_vis_param / l
2965  arg(1) = 2.0_dp * pi * x(1) / l
2966  arg(2) = 2.0_dp * pi * x(2) / l
2967  arg(3) = 2.0_dp * pi * x(3) / l
2968  source = +3.0_dp * ( 2.0_dp * pi / l ) * fact * sin( arg(1) ) * sin( arg(2) ) * sin( arg(3) )
2969  sum = pgm * source
2970  ELSE
2971  sum = 0.0_dp
2972  END IF
2973  ENDIF
2974 
2975  !Calculate RHS VECTOR
2976  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)+sum*rwg
2977  ENDDO !ms
2978  ENDDO !mh
2979  ELSE
2980  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
2981  ENDIF
2982  ENDIF
2983 
2984  ! end: RIGHT HAND SIDE FOR ANALYTIC SOLUTION
2985  !-----------------------------------------------------------------------------------------------------------------------------------
2986 
2987 ! !===================================================================================================================
2988 ! !COUPLING_MATRICES
2989 ! SELECT CASE(EQUATIONS_SET_SUBTYPE)
2990 ! CASE(EQUATIONS_SET_MULTI_COMPARTMENT_DARCY_SUBTYPE,EQUATIONS_SET_ELASTICITY_MULTI_COMPARTMENT_DARCY_INRIA_SUBTYPE)
2991 !
2992 ! !Create FIELD_VARIABLES type, COUPLING_MATRICES type
2993 !
2994 ! !Loop over element rows
2995 ! mhs=0
2996 ! DO mh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
2997 !
2998 ! MESH_COMPONENT_1 = FIELD_VARIABLE%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
2999 ! DEPENDENT_BASIS_1 => DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT_1)%PTR% &
3000 ! & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS
3001 ! QUADRATURE_SCHEME_1 => DEPENDENT_BASIS_1%QUADRATURE% &
3002 ! & QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR
3003 ! RWG = EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(FIELD_U_VARIABLE_TYPE)%PTR%JACOBIAN * &
3004 ! & QUADRATURE_SCHEME_1%GAUSS_WEIGHTS(ng)
3005 !
3006 ! DO ms=1,DEPENDENT_BASIS_1%NUMBER_OF_ELEMENT_PARAMETERS
3007 ! mhs=mhs+1
3008 !
3009 ! DO imatrix=1,Ncompartments
3010 !
3011 ! IF(COUPLING_MATRICES(imatrix)%PTR%UPDATE_MATRIX) THEN
3012 !
3013 ! !Loop over element columns
3014 ! nhs=0
3015 ! ! DO nh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
3016 ! DO nh=1,FIELD_VARIABLES(imatrix)%PTR%NUMBER_OF_COMPONENTS
3017 !
3018 ! MESH_COMPONENT_2 = FIELD_VARIABLE%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
3019 ! DEPENDENT_BASIS_2 => DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT_2)%PTR% &
3020 ! & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS
3021 ! !--- We cannot use two different quadrature schemes here !!!
3022 ! QUADRATURE_SCHEME_2 => DEPENDENT_BASIS_2%QUADRATURE% &
3023 ! & QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR
3024 ! !RWG = EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN * &
3025 ! ! & QUADRATURE_SCHEME_2%GAUSS_WEIGHTS(ng)
3026 !
3027 ! DO ns=1,DEPENDENT_BASIS_2%NUMBER_OF_ELEMENT_PARAMETERS
3028 ! nhs=nhs+1
3029 !
3030 ! !-------------------------------------------------------------------------------------------------------------
3031 ! !velocity test function, velocity trial function
3032 ! !For now, this is only a dummy implementation - this still has to be properly set up.
3033 ! IF(mh==nh.AND.nh<NUMBER_OF_VEL_PRESS_COMPONENTS) THEN
3034 !
3035 ! SUM = 0.0_DP
3036 !
3037 ! PGM=QUADRATURE_SCHEME_1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng)
3038 ! PGN=QUADRATURE_SCHEME_2%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng)
3039 !
3040 ! SUM = SUM + VIS_OVER_PERM_TENSOR( mh, nh ) * PGM * PGN
3041 !
3042 ! COUPLING_MATRICES(imatrix)%PTR%ELEMENT_MATRIX%MATRIX(mhs,nhs) = &
3043 ! & COUPLING_MATRICES(imatrix)%PTR%ELEMENT_MATRIX%MATRIX(mhs,nhs) + SUM * RWG
3044 ! ENDIF
3045 !
3046 ! ENDDO !ns
3047 ! ENDDO !nh
3048 ! ENDIF
3049 ! ENDDO !imatrix
3050 ! ENDDO !ms
3051 ! ENDDO !mh
3052 ! CASE DEFAULT
3053 ! !Do nothing
3054 ! END SELECT
3055  ENDDO !ng
3056 
3057  IF(rhs_vector%UPDATE_VECTOR) THEN
3058  ! Integrate pressure over faces, and add to RHS vector
3059  CALL darcy_finiteelementfaceintegrate(equations_set,element_number,field_variable,err,error,*999)
3060  ENDIF
3061 
3062  ! CHECK STIFFNESS MATRIX WITH CMHEART
3063  IF(diagnostics5) THEN
3064  IF( element_number == 1 ) THEN
3065  ndofs = 0
3066  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3067  mesh_component_1 = field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
3068  dependent_basis_1 => dependent_field%DECOMPOSITION%DOMAIN(mesh_component_1)%PTR% &
3069  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3070  ndofs = ndofs + dependent_basis_1%NUMBER_OF_ELEMENT_PARAMETERS
3071  END DO
3072 
3073  CALL write_string(diagnostic_output_type,"Element Matrix for element number 1 (Darcy):",err,error,*999)
3074  DO mhs=1,ndofs
3075  CALL write_string_value(diagnostic_output_type,"row number = ",mhs,err,error,*999)
3076  CALL write_string_vector(diagnostic_output_type,1,1,ndofs,ndofs,ndofs,&
3077  & stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,:), &
3078  & '("",4(X,E13.6))','4(4(X,E13.6))',err,error,*999)
3079  CALL write_string(diagnostic_output_type," ",err,error,*999)
3080  END DO
3081  END IF
3082  END IF
3083 
3084  !Scale factor adjustment
3085  IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
3086  CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
3087  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
3088  mhs=0
3089  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3090  !Loop over element rows
3091  mesh_component_1=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
3092  dependent_basis_1=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_1)%PTR% &
3093  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3094  DO ms=1,dependent_basis_1%NUMBER_OF_ELEMENT_PARAMETERS
3095  mhs=mhs+1
3096  nhs=0
3097  IF(ASSOCIATED(stiffness_matrix).AND.ASSOCIATED(damping_matrix)) THEN
3098  IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX) THEN
3099  !Loop over element columns
3100  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
3101  mesh_component_2=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
3102  dependent_basis_2=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_2)%PTR% &
3103  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3104  DO ns=1,dependent_basis_2%NUMBER_OF_ELEMENT_PARAMETERS
3105  nhs=nhs+1
3106  IF(stiffness_matrix%UPDATE_MATRIX)THEN
3107  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
3108  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
3109  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
3110  END IF
3111  IF(damping_matrix%UPDATE_MATRIX)THEN
3112  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
3113  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
3114  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
3115  END IF
3116  ENDDO !ns
3117  ENDDO !nh
3118  ENDIF
3119  ENDIF
3120  IF(ASSOCIATED(rhs_vector)) THEN
3121  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
3122  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
3123  ENDIF
3124  IF(ASSOCIATED(source_vector)) THEN
3125  IF(source_vector%UPDATE_VECTOR) source_vector%ELEMENT_VECTOR%VECTOR(mhs)= &
3126  & source_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
3127  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
3128  ENDIF
3129  ENDDO !ms
3130  ENDDO !mh
3131  ENDIF
3132 
3133  ! RESTORE ALL POINTERS CALL PARAMATER_SET_FIELD_DATA_RESTORE
3134 
3135  CASE DEFAULT
3136  local_error="Equations set subtype "//trim(number_to_vstring(equations_set_subtype,"*",err,error))// &
3137  & " is not valid for a Darcy equation type of a fluid mechanics equations set class."
3138  CALL flagerror(local_error,err,error,*999)
3139  END SELECT
3140 
3141  ELSE
3142  CALL flagerror("Equations set equations is not associated.",err,error,*999)
3143  ENDIF
3144  ELSE
3145  CALL flagerror("Equations set is not associated.",err,error,*999)
3146  ENDIF
3147 
3148  exits("DARCY_EQUATION_FINITE_ELEMENT_CALCULATE")
3149  RETURN
3150 999 errorsexits("DARCY_EQUATION_FINITE_ELEMENT_CALCULATE",err,error)
3151  RETURN 1
3152  END SUBROUTINE darcy_equation_finite_element_calculate
3153 
3154  !
3155  !================================================================================================================================
3156  !
3157 
3160  SUBROUTINE darcy_finiteelementfaceintegrate(equationsSet,elementNumber,dependentVariable,err,error,*)
3161 
3162  !Argument variables
3163  TYPE(equations_set_type), POINTER :: equationsSet
3164  INTEGER(INTG), INTENT(IN) :: elementNumber
3165  TYPE(field_variable_type), POINTER :: dependentVariable
3166  INTEGER(INTG), INTENT(OUT) :: err
3167  TYPE(varying_string), INTENT(OUT) :: error
3168  !Local variables
3169  TYPE(decomposition_type), POINTER :: decomposition
3170  TYPE(decomposition_element_type), POINTER :: decompElement
3171  TYPE(basis_type), POINTER :: dependentBasis
3172  TYPE(equations_type), POINTER :: equations
3173  TYPE(equations_matrices_type), POINTER :: equationsMatrices
3174  TYPE(decomposition_face_type), POINTER :: face
3175  TYPE(basis_type), POINTER :: faceBasis
3176  TYPE(field_interpolated_point_type), POINTER :: dependentInterpolatedPoint
3177  TYPE(field_interpolation_parameters_type), POINTER :: dependentInterpolationParameters
3178  TYPE(quadrature_scheme_type), POINTER :: faceQuadratureScheme
3179  TYPE(field_interpolated_point_type), POINTER :: geometricInterpolatedPoint
3180  TYPE(field_interpolation_parameters_type), POINTER :: geometricInterpolationParameters
3181  TYPE(field_interpolated_point_metrics_type), POINTER :: pointMetrics
3182  TYPE(equations_matrices_rhs_type), POINTER :: rhsVector
3183  INTEGER(INTG) :: faceIdx, faceNumber
3184  INTEGER(INTG) :: componentIdx, gaussIdx
3185  INTEGER(INTG) :: elementBaseDofIdx, faceNodeIdx, elementNodeIdx
3186  INTEGER(INTG) :: faceNodeDerivativeIdx, meshComponentNumber, nodeDerivativeIdx, parameterIdx
3187  INTEGER(INTG) :: faceParameterIdx, elementDofIdx, normalComponentIdx
3188  REAL(DP) :: gaussWeight, normalProjection, pressureGauss
3189 
3190  enters("Darcy_FiniteElementFaceIntegrate",err,error,*999)
3191 
3192  NULLIFY(decomposition)
3193  NULLIFY(decompelement)
3194  NULLIFY(dependentbasis)
3195  NULLIFY(equations)
3196  NULLIFY(equationsmatrices)
3197  NULLIFY(face)
3198  NULLIFY(facebasis)
3199  NULLIFY(facequadraturescheme)
3200  NULLIFY(dependentinterpolatedpoint)
3201  NULLIFY(dependentinterpolationparameters)
3202  NULLIFY(geometricinterpolatedpoint)
3203  NULLIFY(geometricinterpolationparameters)
3204  NULLIFY(rhsvector)
3205 
3206  equations=>equationsset%EQUATIONS
3207  IF(ASSOCIATED(equationsset)) THEN
3208  equations=>equationsset%EQUATIONS
3209  IF(ASSOCIATED(equations)) THEN
3210  equationsmatrices=>equations%EQUATIONS_MATRICES
3211  IF(ASSOCIATED(equationsmatrices)) THEN
3212  rhsvector=>equationsmatrices%RHS_VECTOR
3213  END IF
3214  ELSE
3215  CALL flagerror("Equations set equations is not associated.",err,error,*999)
3216  END IF
3217  ELSE
3218  CALL flagerror("Equations set is not associated.",err,error,*999)
3219  END IF
3220 
3221  IF(.NOT.ALLOCATED(equationsset%specification)) THEN
3222  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
3223  ELSE IF(SIZE(equationsset%specification,1)/=3) THEN
3224  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
3225  & err,error,*999)
3226  END IF
3227  SELECT CASE(equationsset%specification(3))
3228  CASE(equations_set_standard_darcy_subtype, &
3229  & equations_set_quasistatic_darcy_subtype, &
3230  & equations_set_ale_darcy_subtype, &
3231  & equations_set_transient_darcy_subtype)
3232 
3233  !Get the mesh decomposition and basis for this element
3234  decomposition=>dependentvariable%FIELD%DECOMPOSITION
3235  !These RHS terms are associated with the equations for the three velocity components,
3236  !rather than the pressure term
3237  meshcomponentnumber=dependentvariable%COMPONENTS(1)%MESH_COMPONENT_NUMBER
3238  dependentbasis=>decomposition%DOMAIN(meshcomponentnumber)%PTR%TOPOLOGY%ELEMENTS% &
3239  & elements(elementnumber)%BASIS
3240  decompelement=>decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(elementnumber)
3241 
3242  !Only add RHS terms if the face geometric parameters are calculated
3243  IF(decomposition%CALCULATE_FACES) THEN
3244  !Get interpolation parameters and point for Darcy pressure
3245  dependentinterpolationparameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS( &
3246  & dependentvariable%VARIABLE_TYPE)%PTR
3247  dependentinterpolatedpoint=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT( &
3248  & dependentvariable%VARIABLE_TYPE)%PTR
3249 
3250  DO faceidx=1,dependentbasis%NUMBER_OF_LOCAL_FACES
3251  !Get the face normal and quadrature information
3252  IF(ALLOCATED(decompelement%ELEMENT_FACES)) THEN
3253  facenumber=decompelement%ELEMENT_FACES(faceidx)
3254  ELSE
3255  CALL flagerror("Decomposition element faces is not allocated.",err,error,*999)
3256  END IF
3257  face=>decomposition%TOPOLOGY%FACES%FACES(facenumber)
3258  !This speeds things up but is also important, as non-boundary faces have an XI_DIRECTION that might
3259  !correspond to the other element.
3260  IF(.NOT.(face%BOUNDARY_FACE)) cycle
3261  CALL field_interpolation_parameters_face_get(field_values_set_type,facenumber, &
3262  & dependentinterpolationparameters,err,error,*999)
3263  normalcomponentidx=abs(face%XI_DIRECTION)
3264  facebasis=>decomposition%DOMAIN(meshcomponentnumber)%PTR%TOPOLOGY%FACES%FACES(facenumber)%BASIS
3265  facequadraturescheme=>facebasis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
3266 
3267  DO gaussidx=1,facequadraturescheme%NUMBER_OF_GAUSS
3268  gaussweight=facequadraturescheme%GAUSS_WEIGHTS(gaussidx)
3269  !Get interpolated Darcy pressure
3270  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,gaussidx, &
3271  & dependentinterpolatedpoint,err,error,*999)
3272  pressuregauss=dependentinterpolatedpoint%values(4,1) !(component,derivative)
3273 
3274  !Use the geometric field to find the face normal and the Jacobian for the face integral
3275  geometricinterpolationparameters=>equations%INTERPOLATION%GEOMETRIC_INTERP_PARAMETERS( &
3276  & field_u_variable_type)%PTR
3277  CALL field_interpolation_parameters_element_get(field_values_set_type,elementnumber, &
3278  & geometricinterpolationparameters,err,error,*999)
3279  geometricinterpolatedpoint=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
3280  CALL field_interpolate_local_face_gauss(first_part_deriv,basis_default_quadrature_scheme,faceidx,gaussidx, &
3281  & geometricinterpolatedpoint,err,error,*999)
3282  !Calculate the metric tensors and Jacobian
3283  pointmetrics=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR
3284  CALL field_interpolated_point_metrics_calculate(coordinate_jacobian_volume_type,pointmetrics,err,error,*999)
3285 
3286  DO componentidx=1,dependentvariable%NUMBER_OF_COMPONENTS-1
3287  normalprojection=dot_product(pointmetrics%GU(normalcomponentidx,:),pointmetrics%DX_DXI(componentidx,:))
3288  IF(face%XI_DIRECTION<0) THEN
3289  normalprojection=-normalprojection
3290  END IF
3291  IF(abs(normalprojection)<zero_tolerance) cycle
3292  !Work out the first index of the rhs vector for this element - 1
3293  elementbasedofidx=dependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(componentidx-1)
3294  DO facenodeidx=1,facebasis%NUMBER_OF_NODES
3295  elementnodeidx=dependentbasis%NODE_NUMBERS_IN_LOCAL_FACE(facenodeidx,faceidx)
3296  DO facenodederivativeidx=1,facebasis%NUMBER_OF_DERIVATIVES(facenodeidx)
3297  nodederivativeidx=dependentbasis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(facenodederivativeidx,facenodeidx,faceidx)
3298  parameteridx=dependentbasis%ELEMENT_PARAMETER_INDEX(nodederivativeidx,elementnodeidx)
3299  faceparameteridx=facebasis%ELEMENT_PARAMETER_INDEX(facenodederivativeidx,facenodeidx)
3300  elementdofidx=elementbasedofidx+parameteridx
3301  rhsvector%ELEMENT_VECTOR%VECTOR(elementdofidx) = rhsvector%ELEMENT_VECTOR%VECTOR(elementdofidx) - &
3302  & gaussweight*pressuregauss*normalprojection* &
3303  & facequadraturescheme%GAUSS_BASIS_FNS(faceparameteridx,no_part_deriv,gaussidx)* &
3304  & pointmetrics%JACOBIAN
3305  END DO !nodeDerivativeIdx
3306  END DO !faceNodeIdx
3307  END DO !componentIdx
3308  END DO !gaussIdx
3309  END DO !faceIdx
3310  END IF !decomposition%calculate_faces
3311 
3312  CASE DEFAULT
3313  ! Do nothing for other equation set subtypes
3314  END SELECT
3315 
3316  exits("Darcy_FiniteElementFaceIntegrate")
3317  RETURN
3318 999 errorsexits("Darcy_FiniteElementFaceIntegrate",err,error)
3319  RETURN 1
3320  END SUBROUTINE
3321 
3322  !
3323  !================================================================================================================================
3324  !
3325 
3327  SUBROUTINE darcy_equationssetspecificationset(equationsSet,specification,err,error,*)
3328 
3329  !Argument variables
3330  TYPE(equations_set_type), POINTER :: equationsSet
3331  INTEGER(INTG), INTENT(IN) :: specification(:)
3332  INTEGER(INTG), INTENT(OUT) :: err
3333  TYPE(varying_string), INTENT(OUT) :: error
3334  !Local Variables
3335  TYPE(varying_string) :: localError
3336  INTEGER(INTG) :: subtype
3337 
3338  enters("Darcy_EquationsSetSpecificationSet",err,error,*999)
3339 
3340  IF(ASSOCIATED(equationsset)) THEN
3341  IF(SIZE(specification,1)/=3) THEN
3342  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
3343  & err,error,*999)
3344  END IF
3345  subtype=specification(3)
3346  SELECT CASE(subtype)
3347  CASE(equations_set_standard_darcy_subtype, &
3348  & equations_set_quasistatic_darcy_subtype, &
3349  & equations_set_ale_darcy_subtype, &
3350  & equations_set_transient_darcy_subtype, &
3351  & equations_set_incompressible_finite_elasticity_darcy_subtype, &
3352  & equations_set_elasticity_darcy_inria_model_subtype, &
3353  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
3354  & equations_set_transient_ale_darcy_subtype, &
3355  & equations_set_multi_compartment_darcy_subtype, &
3356  & equations_set_elasticity_multi_compartment_darcy_inria_subtype, &
3357  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
3358  !ok
3359  CASE DEFAULT
3360  localerror="The third equations set specification of "//trim(numbertovstring(specification(3),"*",err,error))// &
3361  & " is not valid for a Darcy type of a fluid mechanics equations set."
3362  CALL flagerror(localerror,err,error,*999)
3363  END SELECT
3364  !Set full specification
3365  IF(ALLOCATED(equationsset%specification)) THEN
3366  CALL flagerror("Equations set specification is already allocated.",err,error,*999)
3367  ELSE
3368  ALLOCATE(equationsset%specification(3),stat=err)
3369  IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
3370  END IF
3371  equationsset%specification(1:3)=[equations_set_fluid_mechanics_class,equations_set_darcy_equation_type,subtype]
3372  ELSE
3373  CALL flagerror("Equations set is not associated.",err,error,*999)
3374  END IF
3375 
3376  exits("Darcy_EquationsSetSpecificationSet")
3377  RETURN
3378 999 errors("Darcy_EquationsSetSpecificationSet",err,error)
3379  exits("Darcy_EquationsSetSpecificationSet")
3380  RETURN 1
3381 
3382  END SUBROUTINE darcy_equationssetspecificationset
3383 
3384  !
3385  !================================================================================================================================
3386  !
3387 
3389  SUBROUTINE darcy_problemspecificationset(problem,problemSpecification,err,error,*)
3390 
3391  !Argument variables
3392  TYPE(problem_type), POINTER :: problem
3393  INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
3394  INTEGER(INTG), INTENT(OUT) :: err
3395  TYPE(varying_string), INTENT(OUT) :: error
3396  !Local Variables
3397  TYPE(varying_string) :: localError
3398  INTEGER(INTG) :: problemSubtype
3399 
3400  enters("Darcy_ProblemSpecificationSet",err,error,*998)
3401 
3402  IF(ASSOCIATED(problem)) THEN
3403  IF(SIZE(problemspecification,1)==3) THEN
3404  problemsubtype=problemspecification(3)
3405  SELECT CASE(problemsubtype)
3406  CASE(problem_standard_darcy_subtype, &
3407  & problem_quasistatic_darcy_subtype, &
3408  & problem_ale_darcy_subtype, &
3409  & problem_transient_darcy_subtype, &
3410  & problem_pgm_darcy_subtype, &
3411  & problem_pgm_transient_darcy_subtype)
3412  !All ok
3413  CASE DEFAULT
3414  localerror="The third problem subtype of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
3415  & " is not valid for a Darcy type of a fluid mechanics problem."
3416  CALL flagerror(localerror,err,error,*998)
3417  END SELECT
3418  IF(ALLOCATED(problem%specification)) THEN
3419  CALL flagerror("Problem specification is already allocated.",err,error,*998)
3420  ELSE
3421  ALLOCATE(problem%specification(3),stat=err)
3422  IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
3423  END IF
3424  problem%specification(1:3)=[problem_fluid_mechanics_class,problem_darcy_equation_type,problemsubtype]
3425  ELSE
3426  CALL flagerror("Darcy problem specification must have three entries.",err,error,*998)
3427  END IF
3428  ELSE
3429  CALL flagerror("Problem is not associated.",err,error,*998)
3430  END IF
3431 
3432  exits("Darcy_ProblemSpecificationSet")
3433  RETURN
3434 999 IF(ALLOCATED(problem%specification)) DEALLOCATE(problem%specification)
3435 998 errors("Darcy_ProblemSpecificationSet",err,error)
3436  exits("Darcy_ProblemSpecificationSet")
3437  RETURN 1
3438 
3439  END SUBROUTINE darcy_problemspecificationset
3440 
3441  !
3442  !================================================================================================================================
3443  !
3444 
3446 ! SUBROUTINE DARCY_EQUATION_PROBLEM_STANDARD_SETUP(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
3447  SUBROUTINE darcy_equation_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
3448 
3449  !Argument variables
3450  TYPE(problem_type), POINTER :: PROBLEM
3451  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
3452  INTEGER(INTG), INTENT(OUT) :: ERR
3453  TYPE(varying_string), INTENT(OUT) :: ERROR
3454  !Local Variables
3455  TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
3456  TYPE(solver_type), POINTER :: SOLVER, SOLVER_MAT_PROPERTIES
3457  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS, SOLVER_EQUATIONS_MAT_PROPERTIES
3458  TYPE(solvers_type), POINTER :: SOLVERS
3459  TYPE(varying_string) :: LOCAL_ERROR
3460 
3461 ! ENTERS("DARCY_EQUATION_PROBLEM_STANDARD_SETUP",ERR,ERROR,*999)
3462  enters("DARCY_EQUATION_PROBLEM_SETUP",err,error,*999)
3463 
3464  NULLIFY(control_loop)
3465  NULLIFY(solvers)
3466  NULLIFY(solver)
3467  NULLIFY(solver_mat_properties)
3468  NULLIFY(solver_equations)
3469  NULLIFY(solver_equations_mat_properties)
3470  IF(ASSOCIATED(problem)) THEN
3471  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
3472  CALL flagerror("Problem specification is not allocated.",err,error,*999)
3473  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
3474  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
3475  END IF
3476  SELECT CASE(problem%SPECIFICATION(3))
3477 
3478  !-----------------------------------------------------------------
3479  ! s t a n d a r d D a r c y
3480  !-----------------------------------------------------------------
3481  CASE(problem_standard_darcy_subtype)
3482  SELECT CASE(problem_setup%SETUP_TYPE)
3483  CASE(problem_setup_initial_type)
3484  SELECT CASE(problem_setup%ACTION_TYPE)
3485  CASE(problem_setup_start_action)
3486  !Do nothing????
3487  CASE(problem_setup_finish_action)
3488  !Do nothing???
3489  CASE DEFAULT
3490  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3491  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3492  & " is invalid for a standard Darcy equation."
3493  CALL flagerror(local_error,err,error,*999)
3494  END SELECT
3495  CASE(problem_setup_control_type)
3496  SELECT CASE(problem_setup%ACTION_TYPE)
3497  CASE(problem_setup_start_action)
3498  !Set up a simple control loop
3499  CALL control_loop_create_start(problem,control_loop,err,error,*999)
3500  CASE(problem_setup_finish_action)
3501  !Finish the control loops
3502  control_loop_root=>problem%CONTROL_LOOP
3503  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3504  CALL control_loop_create_finish(control_loop,err,error,*999)
3505  CASE DEFAULT
3506  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3507  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3508  & " is invalid for a standard Darcy equation."
3509  CALL flagerror(local_error,err,error,*999)
3510  END SELECT
3511  CASE(problem_setup_solvers_type)
3512  !Get the control loop
3513  control_loop_root=>problem%CONTROL_LOOP
3514  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3515  SELECT CASE(problem_setup%ACTION_TYPE)
3516  CASE(problem_setup_start_action)
3517  !Start the solvers creation
3518  CALL solvers_create_start(control_loop,solvers,err,error,*999)
3519  CALL solvers_number_set(solvers,1,err,error,*999)
3520  !Set the solver to be a linear solver
3521  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3522  CALL solver_type_set(solver,solver_linear_type,err,error,*999)
3523  !Set solver defaults
3524  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
3525  CASE(problem_setup_finish_action)
3526  !Get the solvers
3527  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3528  !Finish the solvers creation
3529  CALL solvers_create_finish(solvers,err,error,*999)
3530  CASE DEFAULT
3531  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3532  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3533  & " is invalid for a standard Darcy equation."
3534  CALL flagerror(local_error,err,error,*999)
3535  END SELECT
3536  CASE(problem_setup_solver_equations_type)
3537  SELECT CASE(problem_setup%ACTION_TYPE)
3538  CASE(problem_setup_start_action)
3539  !Get the control loop
3540  control_loop_root=>problem%CONTROL_LOOP
3541  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3542  !Get the solver
3543  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3544  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3545  !Create the solver equations
3546  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
3547  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
3548  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_static,err,error,*999)
3549  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
3550  CASE(problem_setup_finish_action)
3551  !Get the control loop
3552  control_loop_root=>problem%CONTROL_LOOP
3553  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3554  !Get the solver equations
3555  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3556  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3557  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
3558  !Finish the solver equations creation
3559  CALL solver_equations_create_finish(solver_equations,err,error,*999)
3560  CASE DEFAULT
3561  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3562  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3563  & " is invalid for a standard Darcy equation."
3564  CALL flagerror(local_error,err,error,*999)
3565  END SELECT
3566  CASE DEFAULT
3567  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3568  & " is invalid for a standard Darcy equation."
3569  CALL flagerror(local_error,err,error,*999)
3570  END SELECT
3571 
3572  !-----------------------------------------------------------------
3573  ! q u a s i s t a t i c D a r c y
3574  !-----------------------------------------------------------------
3575  CASE(problem_quasistatic_darcy_subtype)
3576  SELECT CASE(problem_setup%SETUP_TYPE)
3577  CASE(problem_setup_initial_type)
3578  SELECT CASE(problem_setup%ACTION_TYPE)
3579  CASE(problem_setup_start_action)
3580  !Do nothing????
3581  CASE(problem_setup_finish_action)
3582  !Do nothing???
3583  CASE DEFAULT
3584  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3585  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3586  & " is invalid for a quasistatic Darcy equation."
3587  CALL flagerror(local_error,err,error,*999)
3588  END SELECT
3589  CASE(problem_setup_control_type)
3590  SELECT CASE(problem_setup%ACTION_TYPE)
3591  CASE(problem_setup_start_action)
3592  !Set up a time control loop
3593  CALL control_loop_create_start(problem,control_loop,err,error,*999)
3594  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
3595  CASE(problem_setup_finish_action)
3596  !Finish the control loops
3597  control_loop_root=>problem%CONTROL_LOOP
3598  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3599  CALL control_loop_create_finish(control_loop,err,error,*999)
3600  CASE DEFAULT
3601  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3602  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3603  & " is invalid for a quasistatic Darcy equation."
3604  CALL flagerror(local_error,err,error,*999)
3605  END SELECT
3606  CASE(problem_setup_solvers_type)
3607  !Get the control loop
3608  control_loop_root=>problem%CONTROL_LOOP
3609  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3610  SELECT CASE(problem_setup%ACTION_TYPE)
3611  CASE(problem_setup_start_action)
3612  !Start the solvers creation
3613  CALL solvers_create_start(control_loop,solvers,err,error,*999)
3614  CALL solvers_number_set(solvers,1,err,error,*999)
3615  !Set the solver to be a linear solver
3616  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3617  CALL solver_type_set(solver,solver_linear_type,err,error,*999)
3618  !Set solver defaults
3619  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
3620  CASE(problem_setup_finish_action)
3621  !Get the solvers
3622  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3623  !Finish the solvers creation
3624  CALL solvers_create_finish(solvers,err,error,*999)
3625  CASE DEFAULT
3626  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3627  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3628  & " is invalid for a quasistatic Darcy equation."
3629  CALL flagerror(local_error,err,error,*999)
3630  END SELECT
3631  CASE(problem_setup_solver_equations_type)
3632  SELECT CASE(problem_setup%ACTION_TYPE)
3633  CASE(problem_setup_start_action)
3634  !Get the control loop
3635  control_loop_root=>problem%CONTROL_LOOP
3636  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3637  !Get the solver
3638  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3639  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3640  !Create the solver equations
3641  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
3642  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
3643  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_quasistatic,err,error,*999)
3644  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
3645  CASE(problem_setup_finish_action)
3646  !Get the control loop
3647  control_loop_root=>problem%CONTROL_LOOP
3648  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3649  !Get the solver equations
3650  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3651  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3652  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
3653  !Finish the solver equations creation
3654  CALL solver_equations_create_finish(solver_equations,err,error,*999)
3655  CASE DEFAULT
3656  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3657  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3658  & " is invalid for a quasistatic Darcy equation."
3659  CALL flagerror(local_error,err,error,*999)
3660  END SELECT
3661  CASE DEFAULT
3662  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3663  & " is invalid for a quasistatic Darcy equation."
3664  CALL flagerror(local_error,err,error,*999)
3665  END SELECT
3666 
3667  !-----------------------------------------------------------------
3668  ! A L E / P G M D a r c y
3669  !-----------------------------------------------------------------
3670  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype)
3671  SELECT CASE(problem_setup%SETUP_TYPE)
3672  CASE(problem_setup_initial_type)
3673  SELECT CASE(problem_setup%ACTION_TYPE)
3674  CASE(problem_setup_start_action)
3675  !Do nothing????
3676  CASE(problem_setup_finish_action)
3677  !Do nothing???
3678  CASE DEFAULT
3679  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3680  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3681  & " is invalid for an ALE Darcy equation."
3682  CALL flagerror(local_error,err,error,*999)
3683  END SELECT
3684  CASE(problem_setup_control_type)
3685  SELECT CASE(problem_setup%ACTION_TYPE)
3686  CASE(problem_setup_start_action)
3687  !Set up a time control loop
3688  CALL control_loop_create_start(problem,control_loop,err,error,*999)
3689  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
3690  CASE(problem_setup_finish_action)
3691  !Finish the control loops
3692  control_loop_root=>problem%CONTROL_LOOP
3693  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3694  CALL control_loop_create_finish(control_loop,err,error,*999)
3695  CASE DEFAULT
3696  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3697  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3698  & " is invalid for an ALE Darcy equation."
3699  CALL flagerror(local_error,err,error,*999)
3700  END SELECT
3701  CASE(problem_setup_solvers_type)
3702  !Get the control loop
3703  control_loop_root=>problem%CONTROL_LOOP
3704  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3705  SELECT CASE(problem_setup%ACTION_TYPE)
3706  CASE(problem_setup_start_action)
3707  !Start the solvers creation
3708  CALL solvers_create_start(control_loop,solvers,err,error,*999)
3709  CALL solvers_number_set(solvers,2,err,error,*999)
3710  !
3711  !Set the first solver to be a linear solver for the material update
3712  CALL solvers_solver_get(solvers,1,solver_mat_properties,err,error,*999)
3713  CALL solver_type_set(solver_mat_properties,solver_linear_type,err,error,*999)
3714  CALL solver_library_type_set(solver_mat_properties,solver_petsc_library,err,error,*999)
3715  !
3716  !Set the second solver to be a linear solver for the ALE Darcy
3717  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
3718  CALL solver_type_set(solver,solver_linear_type,err,error,*999)
3719  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
3720  CASE(problem_setup_finish_action)
3721  !Get the solvers
3722  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3723  !Finish the solvers creation
3724  CALL solvers_create_finish(solvers,err,error,*999)
3725  CASE DEFAULT
3726  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3727  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3728  & " is invalid for an ALE Darcy equation."
3729  CALL flagerror(local_error,err,error,*999)
3730  END SELECT
3731  CASE(problem_setup_solver_equations_type)
3732  SELECT CASE(problem_setup%ACTION_TYPE)
3733  CASE(problem_setup_start_action)
3734  !Get the control loop and solvers
3735  control_loop_root=>problem%CONTROL_LOOP
3736  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3737  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3738  !Get the material-properties solver and create the material-properties solver equations
3739  CALL solvers_solver_get(solvers,1,solver_mat_properties,err,error,*999)
3740  CALL solver_equations_create_start(solver_mat_properties,solver_equations_mat_properties,err,error,*999)
3741  CALL solver_equations_linearity_type_set(solver_equations_mat_properties,solver_equations_linear,err,error,*999)
3742  CALL solver_equations_time_dependence_type_set(solver_equations_mat_properties,solver_equations_quasistatic, &
3743  & err,error,*999)
3744  CALL solver_equations_sparsity_type_set(solver_equations_mat_properties,solver_sparse_matrices,err,error,*999)
3745  !Get the Darcy-ALE solver and create the Darcy-ALE solver equations
3746  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
3747  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
3748  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
3749  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_quasistatic,err,error,*999)
3750  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
3751  CASE(problem_setup_finish_action)
3752  !Get the control loop
3753  control_loop_root=>problem%CONTROL_LOOP
3754  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3755  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3756  !Finish the creation of the material-properties solver equations
3757  CALL solvers_solver_get(solvers,1,solver_mat_properties,err,error,*999)
3758  CALL solver_solver_equations_get(solver_mat_properties,solver_equations_mat_properties,err,error,*999)
3759  CALL solver_equations_create_finish(solver_equations_mat_properties,err,error,*999)
3760  !Finish the creation of the Darcy-ALE solver equations
3761  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
3762  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
3763  CALL solver_equations_create_finish(solver_equations,err,error,*999)
3764  CASE DEFAULT
3765  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3766  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3767  & " is invalid for an ALE Darcy equation."
3768  CALL flagerror(local_error,err,error,*999)
3769  END SELECT
3770  CASE DEFAULT
3771  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3772  & " is invalid for an ALE Darcy equation."
3773  CALL flagerror(local_error,err,error,*999)
3774  END SELECT
3775 
3776  !-----------------------------------------------------------------
3777  ! D Y N A M I C A L E / P G M D a r c y
3778  !-----------------------------------------------------------------
3779  CASE(problem_pgm_transient_darcy_subtype)
3780  SELECT CASE(problem_setup%SETUP_TYPE)
3781  CASE(problem_setup_initial_type)
3782  SELECT CASE(problem_setup%ACTION_TYPE)
3783  CASE(problem_setup_start_action)
3784  !Do nothing????
3785  CASE(problem_setup_finish_action)
3786  !Do nothing???
3787  CASE DEFAULT
3788  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3789  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3790  & " is invalid for an ALE Darcy equation."
3791  CALL flagerror(local_error,err,error,*999)
3792  END SELECT
3793  CASE(problem_setup_control_type)
3794  SELECT CASE(problem_setup%ACTION_TYPE)
3795  CASE(problem_setup_start_action)
3796  !Set up a time control loop
3797  CALL control_loop_create_start(problem,control_loop,err,error,*999)
3798  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
3799  CASE(problem_setup_finish_action)
3800  !Finish the control loops
3801  control_loop_root=>problem%CONTROL_LOOP
3802  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3803  CALL control_loop_create_finish(control_loop,err,error,*999)
3804  CASE DEFAULT
3805  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3806  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3807  & " is invalid for an ALE Darcy equation."
3808  CALL flagerror(local_error,err,error,*999)
3809  END SELECT
3810  CASE(problem_setup_solvers_type)
3811  !Get the control loop
3812  control_loop_root=>problem%CONTROL_LOOP
3813  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3814  SELECT CASE(problem_setup%ACTION_TYPE)
3815  CASE(problem_setup_start_action)
3816  !Start the solvers creation
3817  CALL solvers_create_start(control_loop,solvers,err,error,*999)
3818  CALL solvers_number_set(solvers,2,err,error,*999)
3819  !
3820  !Set the first solver to be a linear solver for the material update
3821  CALL solvers_solver_get(solvers,1,solver_mat_properties,err,error,*999)
3822  CALL solver_type_set(solver_mat_properties,solver_linear_type,err,error,*999)
3823  CALL solver_library_type_set(solver_mat_properties,solver_petsc_library,err,error,*999)
3824  !
3825  !Set the second solver to be a first order dynamic solver for the ALE Darcy
3826  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
3827  CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
3828  CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
3829  !Set solver defaults
3830  CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
3831  CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
3832  CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
3833 ! CALL SOLVER_LIBRARY_TYPE_SET(SOLVER,SOLVER_PETSC_LIBRARY,ERR,ERROR,*999)
3834  CASE(problem_setup_finish_action)
3835  !Get the solvers
3836  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3837  !Finish the solvers creation
3838  CALL solvers_create_finish(solvers,err,error,*999)
3839  CASE DEFAULT
3840  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3841  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3842  & " is invalid for an ALE Darcy equation."
3843  CALL flagerror(local_error,err,error,*999)
3844  END SELECT
3845  CASE(problem_setup_solver_equations_type)
3846  SELECT CASE(problem_setup%ACTION_TYPE)
3847  CASE(problem_setup_start_action)
3848  !Get the control loop and solvers
3849  control_loop_root=>problem%CONTROL_LOOP
3850  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3851  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3852  !Get the material-properties solver and create the material-properties solver equations
3853  CALL solvers_solver_get(solvers,1,solver_mat_properties,err,error,*999)
3854  CALL solver_equations_create_start(solver_mat_properties,solver_equations_mat_properties,err,error,*999)
3855  CALL solver_equations_linearity_type_set(solver_equations_mat_properties,solver_equations_linear,err,error,*999)
3856  CALL solver_equations_time_dependence_type_set(solver_equations_mat_properties,solver_equations_quasistatic, &
3857  & err,error,*999)
3858  CALL solver_equations_sparsity_type_set(solver_equations_mat_properties,solver_sparse_matrices,err,error,*999)
3859  !Get the Darcy-ALE solver and create the Darcy-ALE solver equations
3860  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
3861  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
3862  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
3863  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,err,error,*999)
3864  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
3865  CASE(problem_setup_finish_action)
3866  !Get the control loop
3867  control_loop_root=>problem%CONTROL_LOOP
3868  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3869  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3870  !Finish the creation of the material-properties solver equations
3871  CALL solvers_solver_get(solvers,1,solver_mat_properties,err,error,*999)
3872  CALL solver_solver_equations_get(solver_mat_properties,solver_equations_mat_properties,err,error,*999)
3873  CALL solver_equations_create_finish(solver_equations_mat_properties,err,error,*999)
3874  !Finish the creation of the Darcy-ALE solver equations
3875  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
3876  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
3877  CALL solver_equations_create_finish(solver_equations,err,error,*999)
3878  CASE DEFAULT
3879  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3880  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3881  & " is invalid for an ALE Darcy equation."
3882  CALL flagerror(local_error,err,error,*999)
3883  END SELECT
3884  CASE DEFAULT
3885  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3886  & " is invalid for an ALE Darcy equation."
3887  CALL flagerror(local_error,err,error,*999)
3888  END SELECT
3889 
3890  !-----------------------------------------------------------------
3891  ! d y n a m i c D a r c y
3892  !-----------------------------------------------------------------
3893  CASE(problem_transient_darcy_subtype)
3894  SELECT CASE(problem_setup%SETUP_TYPE)
3895  CASE(problem_setup_initial_type)
3896  SELECT CASE(problem_setup%ACTION_TYPE)
3897  CASE(problem_setup_start_action)
3898  !Do nothing????
3899  CASE(problem_setup_finish_action)
3900  !Do nothing????
3901  CASE DEFAULT
3902  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3903  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3904  & " is invalid for a transient Darcy fluid."
3905  CALL flagerror(local_error,err,error,*999)
3906  END SELECT
3907  CASE(problem_setup_control_type)
3908  SELECT CASE(problem_setup%ACTION_TYPE)
3909  CASE(problem_setup_start_action)
3910  !Set up a time control loop
3911  CALL control_loop_create_start(problem,control_loop,err,error,*999)
3912  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
3913  CASE(problem_setup_finish_action)
3914  !Finish the control loops
3915  control_loop_root=>problem%CONTROL_LOOP
3916  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3917  CALL control_loop_create_finish(control_loop,err,error,*999)
3918  CASE DEFAULT
3919  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3920  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3921  & " is invalid for a transient Darcy fluid."
3922  CALL flagerror(local_error,err,error,*999)
3923  END SELECT
3924  CASE(problem_setup_solvers_type)
3925  !Get the control loop
3926  control_loop_root=>problem%CONTROL_LOOP
3927  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3928  SELECT CASE(problem_setup%ACTION_TYPE)
3929  CASE(problem_setup_start_action)
3930  !Start the solvers creation
3931  CALL solvers_create_start(control_loop,solvers,err,error,*999)
3932  CALL solvers_number_set(solvers,1,err,error,*999)
3933  !Set the solver to be a first order dynamic solver
3934  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3935  CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
3936  CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
3937  !Set solver defaults
3938  CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
3939  CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
3940  CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
3941  CASE(problem_setup_finish_action)
3942  !Get the solvers
3943  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3944  !Finish the solvers creation
3945  CALL solvers_create_finish(solvers,err,error,*999)
3946  CASE DEFAULT
3947  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3948  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3949  & " is invalid for a transient Darcy fluid."
3950  CALL flagerror(local_error,err,error,*999)
3951  END SELECT
3952  CASE(problem_setup_solver_equations_type)
3953  SELECT CASE(problem_setup%ACTION_TYPE)
3954  CASE(problem_setup_start_action)
3955  !Get the control loop
3956  control_loop_root=>problem%CONTROL_LOOP
3957  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3958  !Get the solver
3959  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3960  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3961  !Create the solver equations
3962  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
3963  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
3964  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,&
3965  & err,error,*999)
3966  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
3967  CASE(problem_setup_finish_action)
3968  !Get the control loop
3969  control_loop_root=>problem%CONTROL_LOOP
3970  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
3971  !Get the solver equations
3972  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
3973  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
3974  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
3975  !Finish the solver equations creation
3976  CALL solver_equations_create_finish(solver_equations,err,error,*999)
3977  CASE DEFAULT
3978  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
3979  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3980  & " is invalid for a transient Darcy fluid."
3981  CALL flagerror(local_error,err,error,*999)
3982  END SELECT
3983  CASE DEFAULT
3984  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
3985  & " is invalid for a transient Darcy fluid."
3986  CALL flagerror(local_error,err,error,*999)
3987  END SELECT
3988 
3989  !-----------------------------------------------------------------
3990  ! c a s e d e f a u l t
3991  !-----------------------------------------------------------------
3992  CASE DEFAULT
3993  local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
3994  & " does not equal a standard, quasistatic or ALE Darcy equation subtype."
3995  CALL flagerror(local_error,err,error,*999)
3996 
3997  END SELECT
3998  ELSE
3999  CALL flagerror("Problem is not associated.",err,error,*999)
4000  ENDIF
4001 
4002 ! EXITS("DARCY_EQUATION_PROBLEM_STANDARD_SETUP")
4003  exits("DARCY_EQUATION_PROBLEM_SETUP")
4004  RETURN
4005 ! 999 ERRORSEXITS("DARCY_EQUATION_PROBLEM_STANDARD_SETUP",ERR,ERROR)
4006 999 errorsexits("DARCY_EQUATION_PROBLEM_SETUP",err,error)
4007  RETURN 1
4008 ! END SUBROUTINE DARCY_EQUATION_PROBLEM_STANDARD_SETUP
4009  END SUBROUTINE darcy_equation_problem_setup
4010 
4011  !
4012  !================================================================================================================================
4013  !
4014 
4016  SUBROUTINE darcy_equation_pre_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4017 
4018  !Argument variables
4019  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4020  TYPE(solver_type), POINTER :: SOLVER
4021  INTEGER(INTG), INTENT(OUT) :: ERR
4022  TYPE(varying_string), INTENT(OUT) :: ERROR
4023 
4024  !Local Variables
4025  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4026  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
4027  TYPE(solver_matrices_type), POINTER :: SOLVER_MATRICES
4028  TYPE(solver_matrix_type), POINTER :: SOLVER_MATRIX
4029  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4030  TYPE(equations_type), POINTER :: EQUATIONS
4031 
4032  TYPE(solver_type), POINTER :: SOLVER_ALE_DARCY
4033  TYPE(varying_string) :: LOCAL_ERROR
4034 
4035  INTEGER(INTG) :: solver_matrix_idx
4036 
4037 
4038  enters("DARCY_EQUATION_PRE_SOLVE",err,error,*999)
4039 
4040  IF(ASSOCIATED(control_loop)) THEN
4041  IF(ASSOCIATED(solver)) THEN
4042  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4043  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4044  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4045  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4046  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
4047  END IF
4048  !--- Set 'SOLVER_NUMBER' depending on CONTROL_LOOP%PROBLEM%SPECIFICATION(3)
4049  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4050  CASE(problem_standard_darcy_subtype,problem_quasistatic_darcy_subtype,problem_transient_darcy_subtype)
4051  solver_number_darcy=1
4052  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype,problem_pgm_transient_darcy_subtype)
4053  solver_number_mat_properties=1
4054  solver_number_darcy=2
4055  CASE(problem_quasistatic_elasticity_transient_darcy_subtype)
4056  solver_number_solid=1
4057  solver_number_darcy=1
4058  CASE(problem_standard_elasticity_darcy_subtype,problem_pgm_elasticity_darcy_subtype, &
4059  & problem_quasistatic_elast_trans_darcy_mat_solve_subtype)
4060  solver_number_solid=1
4061  solver_number_mat_properties=1
4062  solver_number_darcy=2
4063  END SELECT
4064 
4065  !--- Set explicitly 'SOLVER_MATRIX%UPDATE_MATRIX=.TRUE.'
4066  solver_equations=>solver%SOLVER_EQUATIONS
4067  IF(ASSOCIATED(solver_equations)) THEN
4068  solver_mapping=>solver_equations%SOLVER_MAPPING
4069  IF(ASSOCIATED(solver_mapping)) THEN
4070  solver_matrices=>solver_equations%SOLVER_MATRICES
4071  IF(ASSOCIATED(solver_matrices)) THEN
4072  DO solver_matrix_idx=1,solver_mapping%NUMBER_OF_SOLVER_MATRICES
4073  solver_matrix=>solver_matrices%MATRICES(solver_matrix_idx)%PTR
4074  IF(ASSOCIATED(solver_matrix)) THEN
4075  solver_matrix%UPDATE_MATRIX=.true.
4076  ELSE
4077  CALL flagerror("Solver Matrix is not associated.",err,error,*999)
4078  ENDIF
4079  ENDDO
4080  ELSE
4081  CALL flagerror("Solver Matrices is not associated.",err,error,*999)
4082  ENDIF
4083  ELSE
4084  CALL flagerror("Solver mapping is not associated.",err,error,*999)
4085  ENDIF
4086  ELSE
4087  CALL flagerror("Solver equations is not associated.",err,error,*999)
4088  ENDIF
4089 
4090  !--- pre_solve calls for various actions
4091  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4092  CASE(problem_standard_darcy_subtype)
4093  ! do nothing
4094  CASE(problem_quasistatic_darcy_subtype)
4095  ! do nothing
4096  CASE(problem_transient_darcy_subtype)
4097  CALL darcy_presolveupdateboundaryconditions(control_loop,solver,err,error,*999)
4098  CASE(problem_quasistatic_elasticity_transient_darcy_subtype)
4099 
4100  IF((control_loop%LOOP_TYPE==problem_control_simple_type.OR.control_loop%LOOP_TYPE==problem_control_time_loop_type) &
4101  & .AND.solver%GLOBAL_NUMBER==solver_number_darcy) THEN
4102  !--- flags to ensure once-per-time-step output in conjunction with diagnostics
4103  idebug1 = .true.
4104  idebug2 = .true.
4105  idebug3 = .true.
4106 
4107  NULLIFY(solver_ale_darcy)
4108  CALL solvers_solver_get(solver%SOLVERS,solver_number_darcy,solver_ale_darcy,err,error,*999)
4109  equations=>solver_ale_darcy%SOLVER_EQUATIONS%SOLVER_MAPPING%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
4110  IF(ASSOCIATED(equations)) THEN
4111  equations_set=>equations%EQUATIONS_SET
4112  IF(ASSOCIATED(equations_set)) THEN
4113  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
4114  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_incomp_elast_darcy_analytic_darcy)THEN
4115  !call only analytic update and DO NOT call the other standard pre-solve routines as the mesh does not require deformation
4116  CALL darcy_presolveupdateanalyticvalues(control_loop,solver_ale_darcy,err,error,*999)
4117  ENDIF
4118  ELSE
4119  !default
4120  !--- 1.1 Transfer solid displacement to Darcy geometric field
4121  CALL darcy_presolvegetsoliddisplacement(control_loop,solver_ale_darcy,err,error,*999)
4122 
4123  !--- 1.2 Update the mesh (and calculate boundary velocities) PRIOR to solving for new material properties
4124  CALL darcy_equation_pre_solve_ale_update_mesh(control_loop,solver_ale_darcy,err,error,*999)
4125 
4126  ! ! ! !i n p r i n c i p l e c u r r e n t l y d o n o t n e e d t o u p d a t e B C s
4127  ! ! ! !unless:
4128  ! ! ! !--- 1.3 Apply both normal and moving mesh boundary conditions, OR:
4129  ! ! ! !--- 1.3 (Iteratively) Render the boundary impermeable (ellipsoid, general curvilinear mesh)
4130  ! ! ! CALL Darcy_PreSolveUpdateBoundaryConditions(CONTROL_LOOP,SOLVER_ALE_DARCY,ERR,ERROR,*999)
4131  ENDIF
4132  ENDIF
4133  ENDIF
4134  END IF
4135  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype,problem_standard_elasticity_darcy_subtype, &
4136  & problem_pgm_elasticity_darcy_subtype,problem_pgm_transient_darcy_subtype, &
4137  & problem_quasistatic_elast_trans_darcy_mat_solve_subtype)
4138 
4139  IF((control_loop%LOOP_TYPE==problem_control_simple_type.OR.control_loop%LOOP_TYPE==problem_control_time_loop_type) &
4140  & .AND.solver%GLOBAL_NUMBER==solver_number_mat_properties) THEN
4141  !--- flags to ensure once-per-time-step output in conjunction with diagnostics
4142  idebug1 = .true.
4143  idebug2 = .true.
4144  idebug3 = .true.
4145 
4146  NULLIFY(solver_ale_darcy)
4147  CALL solvers_solver_get(solver%SOLVERS,solver_number_darcy,solver_ale_darcy,err,error,*999)
4148  equations=>solver_ale_darcy%SOLVER_EQUATIONS%SOLVER_MAPPING%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
4149  IF(ASSOCIATED(equations)) THEN
4150  equations_set=>equations%EQUATIONS_SET
4151  IF(ASSOCIATED(equations_set)) THEN
4152  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
4153  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_incomp_elast_darcy_analytic_darcy)THEN
4154  !call only analytic update and DO NOT call the other standard pre-solve routines as the mesh does not require deformation
4155  CALL darcy_presolveupdateanalyticvalues(control_loop,solver_ale_darcy,err,error,*999)
4156  ENDIF
4157  ELSE
4158  !default
4159  !--- 1.1 Transfer solid displacement to Darcy geometric field
4160  CALL darcy_presolvegetsoliddisplacement(control_loop,solver_ale_darcy,err,error,*999)
4161 
4162  !--- 1.2 Update the mesh (and calculate boundary velocities) PRIOR to solving for new material properties
4163  CALL darcy_equation_pre_solve_ale_update_mesh(control_loop,solver_ale_darcy,err,error,*999)
4164 
4165  ! ! ! !i n p r i n c i p l e c u r r e n t l y d o n o t n e e d t o u p d a t e B C s
4166  ! ! ! !--- 1.3 Apply both normal and moving mesh boundary conditions
4167  ! ! ! CALL Darcy_PreSolveUpdateBoundaryConditions(CONTROL_LOOP,SOLVER_ALE_DARCY,ERR,ERROR,*999)
4168  ENDIF
4169  ENDIF
4170  ENDIF
4171  ELSE IF((control_loop%LOOP_TYPE==problem_control_simple_type.OR. &
4172  & control_loop%LOOP_TYPE==problem_control_time_loop_type).AND.solver%GLOBAL_NUMBER==solver_number_darcy) THEN
4173 ! ! ! !n o t f o r n o w ! ! !
4174 ! ! ! !--- 2.1 Update the material field
4175 ! ! ! CALL Darcy_PreSolveUpdateMatrixProperties(CONTROL_LOOP,SOLVER,ERR,ERROR,*999)
4176  END IF
4177  CASE DEFAULT
4178  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4179  & " is not valid for a Darcy fluid type of a fluid mechanics problem class."
4180  CALL flagerror(local_error,err,error,*999)
4181  END SELECT
4182  ELSE
4183  CALL flagerror("Problem is not associated.",err,error,*999)
4184  ENDIF
4185  ELSE
4186  CALL flagerror("Solver is not associated.",err,error,*999)
4187  ENDIF
4188  ELSE
4189  CALL flagerror("Control loop is not associated.",err,error,*999)
4190  ENDIF
4191 
4192  exits("DARCY_EQUATION_PRE_SOLVE")
4193  RETURN
4194 999 errorsexits("DARCY_EQUATION_PRE_SOLVE",err,error)
4195  RETURN 1
4196  END SUBROUTINE darcy_equation_pre_solve
4197 
4198  !
4199  !================================================================================================================================
4200  !
4201 
4202  SUBROUTINE darcy_control_time_loop_pre_loop(CONTROL_LOOP,ERR,ERROR,*)
4203 
4204  !Argument variables
4205  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4206  INTEGER(INTG), INTENT(OUT) :: ERR
4207  TYPE(varying_string), INTENT(OUT) :: ERROR
4208 
4209  !Local Variables
4210  TYPE(control_loop_type), POINTER :: CONTROL_LOOP_DARCY
4211  TYPE(solver_type), POINTER :: SOLVER_DARCY
4212 
4213  enters("DARCY_CONTROL_TIME_LOOP_PRE_LOOP",err,error,*999)
4214 
4215  !Get the solver for the Darcy problem
4216  NULLIFY(solver_darcy)
4217  NULLIFY(control_loop_darcy)
4218  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4219  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4220  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4221  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
4222  END IF
4223  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4224  CASE(problem_standard_darcy_subtype,problem_quasistatic_darcy_subtype,problem_transient_darcy_subtype)
4225  !SOLVER_NUMBER_DARCY has to be set here so that store_reference_data and store_previous_data have access to it
4226  solver_number_darcy=1
4227  CALL control_loop_get(control_loop,[control_loop_node],control_loop_darcy,err,error,*999)
4228  CALL solvers_solver_get(control_loop_darcy%SOLVERS,solver_number_darcy,solver_darcy,err,error,*999)
4229  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype,problem_pgm_transient_darcy_subtype)
4230  solver_number_mat_properties=1
4231  solver_number_darcy=2
4232  CALL control_loop_get(control_loop,[control_loop_node],control_loop_darcy,err,error,*999)
4233  CALL solvers_solver_get(control_loop_darcy%SOLVERS,solver_number_darcy,solver_darcy,err,error,*999)
4234  CASE(problem_standard_elasticity_darcy_subtype,problem_pgm_elasticity_darcy_subtype)
4235  solver_number_solid=1
4236  solver_number_mat_properties=1
4237  solver_number_darcy=2
4238  CALL control_loop_get(control_loop,[2,control_loop_node],control_loop_darcy,err,error,*999)
4239  CALL solvers_solver_get(control_loop_darcy%SOLVERS,solver_number_darcy,solver_darcy,err,error,*999)
4240  CASE(problem_quasistatic_elasticity_transient_darcy_subtype)
4241  solver_number_solid=1
4242  solver_number_darcy=1
4243  CALL control_loop_get(control_loop,[1,2,control_loop_node],control_loop_darcy,err,error,*999)
4244  CALL solvers_solver_get(control_loop_darcy%SOLVERS,solver_number_darcy,solver_darcy,err,error,*999)
4245  CASE(problem_quasistatic_elast_trans_darcy_mat_solve_subtype)
4246  solver_number_solid=1
4247  solver_number_mat_properties=1
4248  solver_number_darcy=2
4249  CALL control_loop_get(control_loop,[1,2,control_loop_node],control_loop_darcy,err,error,*999)
4250  CALL solvers_solver_get(control_loop_darcy%SOLVERS,solver_number_darcy,solver_darcy,err,error,*999)
4251  END SELECT
4252 
4253  !If this is the first time step then store reference data
4254  IF(control_loop%TIME_LOOP%ITERATION_NUMBER==0) THEN
4255  IF(control_loop%OUTPUT_TYPE>=control_loop_progress_output) THEN
4256  CALL write_string(general_output_type,'== Storing reference data',err,error,*999)
4257  ENDIF
4258  CALL darcy_presolvestorereferencedata(control_loop,solver_darcy,err,error,*999)
4259  ENDIF
4260 
4261  !Store data of previous time step (mesh position); executed once per time step before subiteration
4262  IF(control_loop%OUTPUT_TYPE>=control_loop_progress_output) THEN
4263  CALL write_string(general_output_type,'== Storing previous data',err,error,*999)
4264  ENDIF
4265  CALL darcy_presolvestorepreviousdata(control_loop,solver_darcy,err,error,*999)
4266 
4267  exits("DARCY_CONTROL_TIME_LOOP_PRE_LOOP")
4268  RETURN
4269 999 errorsexits("DARCY_CONTROL_TIME_LOOP_PRE_LOOP",err,error)
4270  RETURN 1
4271  END SUBROUTINE darcy_control_time_loop_pre_loop
4272 
4273  !
4274  !================================================================================================================================
4275  !
4276 
4278  SUBROUTINE darcy_presolvestorereferencedata(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4279 
4280  !Argument variables
4281  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4282  TYPE(solver_type), POINTER :: SOLVER
4283  INTEGER(INTG), INTENT(OUT) :: ERR
4284  TYPE(varying_string), INTENT(OUT) :: ERROR
4285  !Local Variables
4286  TYPE(field_type), POINTER :: DEPENDENT_FIELD, GEOMETRIC_FIELD
4287  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4288  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
4289  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4290  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
4291  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
4292  TYPE(varying_string) :: LOCAL_ERROR
4293 
4294  REAL(DP) :: ALPHA
4295  REAL(DP), POINTER :: INITIAL_VALUES(:)
4296 
4297  INTEGER(INTG) :: FIELD_VAR_TYPE
4298  INTEGER(INTG) :: NDOFS_TO_PRINT,equations_set_idx
4299 
4300 
4301  enters("Darcy_PreSolveStoreReferenceData",err,error,*999)
4302 
4303  NULLIFY(solver_equations)
4304  NULLIFY(solver_mapping)
4305  NULLIFY(equations_set)
4306 
4307  IF(ASSOCIATED(control_loop)) THEN
4308  IF(ASSOCIATED(solver)) THEN
4309  IF(solver%GLOBAL_NUMBER==solver_number_darcy) THEN
4310  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4311  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4312  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4313  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4314  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
4315  END IF
4316  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4317  CASE(problem_standard_darcy_subtype)
4318  ! do nothing
4319  CASE(problem_quasistatic_darcy_subtype)
4320  ! do nothing
4321  CASE(problem_transient_darcy_subtype)
4322  ! do nothing
4323  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype,problem_standard_elasticity_darcy_subtype, &
4324  & problem_pgm_elasticity_darcy_subtype,problem_pgm_transient_darcy_subtype, &
4325  & problem_quasistatic_elasticity_transient_darcy_subtype,problem_quasistatic_elast_trans_darcy_mat_solve_subtype)
4326  solver_equations=>solver%SOLVER_EQUATIONS
4327  IF(ASSOCIATED(solver_equations)) THEN
4328  solver_mapping=>solver_equations%SOLVER_MAPPING
4329  IF(ASSOCIATED(solver_mapping)) THEN
4330  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
4331  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
4332  IF(ASSOCIATED(equations_set)) THEN
4333  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
4334  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
4335  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
4336  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
4337  & err,error,*999)
4338  END IF
4339  SELECT CASE(equations_set%SPECIFICATION(3))
4340  CASE(equations_set_standard_darcy_subtype)
4341  ! do nothing
4342  CASE(equations_set_quasistatic_darcy_subtype)
4343  ! do nothing
4344  CASE(equations_set_transient_darcy_subtype)
4345  ! do nothing
4346  CASE(equations_set_ale_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
4347  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
4348  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
4349  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
4350  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
4351  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
4352  IF(ASSOCIATED(dependent_field).AND.ASSOCIATED(geometric_field)) THEN
4353  !--- Store the initial (= reference) GEOMETRY field values
4354  alpha = 1.0_dp
4355  CALL field_parameter_sets_copy(geometric_field,field_u_variable_type, &
4356  & field_values_set_type,field_initial_values_set_type,alpha,err,error,*999)
4357  equations_mapping=>equations_set%EQUATIONS%EQUATIONS_MAPPING
4358  IF(ASSOCIATED(equations_mapping)) THEN
4359 
4360  SELECT CASE(equations_set%SPECIFICATION(3))
4361  CASE(equations_set_ale_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype)
4362  field_variable=>equations_mapping%LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
4363  ! '1' associated with linear matrix
4364  CASE(equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
4365  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
4366  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
4367  field_variable=>equations_mapping%DYNAMIC_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
4368  END SELECT
4369 
4370  IF(ASSOCIATED(field_variable)) THEN
4371  field_var_type=field_variable%VARIABLE_TYPE
4372  !--- Store the initial DEPENDENT field values
4373  alpha = 1.0_dp
4374  CALL field_parameter_sets_copy(dependent_field,field_var_type, &
4375  & field_values_set_type,field_initial_values_set_type,alpha,err,error,*999)
4376 
4377  IF(diagnostics1) THEN
4378  NULLIFY(initial_values)
4379  CALL field_parameter_set_data_get(dependent_field,field_var_type, &
4380  & field_initial_values_set_type,initial_values,err,error,*999)
4381  ndofs_to_print = SIZE(initial_values,1)
4382  CALL write_string_vector(diagnostic_output_type,1,1,ndofs_to_print,ndofs_to_print,ndofs_to_print,&
4383  & initial_values, &
4384  & '(" DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE,FIELD_INITIAL_VALUES_SET_TYPE = ",4(X,E13.6))', &
4385  & '4(4(X,E13.6))',err,error,*999)
4386  CALL field_parameter_set_data_restore(dependent_field,field_var_type, &
4387  & field_initial_values_set_type,initial_values,err,error,*999)
4388  ENDIF
4389  ELSE
4390  CALL flagerror("FIELD_VAR_TYPE is not associated.",err,error,*999)
4391  ENDIF
4392  ELSE
4393  CALL flagerror("EQUATIONS_MAPPING is not associated.",err,error,*999)
4394  ENDIF
4395  ELSE
4396  CALL flagerror("Dependent field and / or geometric field is / are not associated.",err,error,*999)
4397  ENDIF
4398  CASE DEFAULT
4399  local_error="Equations set subtype " &
4400  & //trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
4401  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4402  CALL flagerror(local_error,err,error,*999)
4403  END SELECT
4404  ELSE
4405  CALL flagerror("Equations set is not associated.",err,error,*999)
4406  ENDIF
4407  ENDDO
4408  ELSE
4409  CALL flagerror("Solver mapping is not associated.",err,error,*999)
4410  ENDIF
4411  ELSE
4412  CALL flagerror("Solver equations is not associated.",err,error,*999)
4413  ENDIF
4414  CASE DEFAULT
4415  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4416  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4417  CALL flagerror(local_error,err,error,*999)
4418  END SELECT
4419  ELSE
4420  CALL flagerror("Problem is not associated.",err,error,*999)
4421  ENDIF
4422  ELSE
4423  ! do nothing
4424  ENDIF
4425  ELSE
4426  CALL flagerror("Solver is not associated.",err,error,*999)
4427  ENDIF
4428  ELSE
4429  CALL flagerror("Control loop is not associated.",err,error,*999)
4430  ENDIF
4431 
4432 
4433  exits("Darcy_PreSolveStoreReferenceData")
4434  RETURN
4435 999 errorsexits("Darcy_PreSolveStoreReferenceData",err,error)
4436  RETURN 1
4437 
4438  END SUBROUTINE darcy_presolvestorereferencedata
4439 
4440  !
4441  !================================================================================================================================
4442  !
4443 
4445  SUBROUTINE darcy_presolvestorepreviousdata(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4446 
4447  !Argument variables
4448  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4449  TYPE(solver_type), POINTER :: SOLVER
4450  INTEGER(INTG), INTENT(OUT) :: ERR
4451  TYPE(varying_string), INTENT(OUT) :: ERROR
4452  !Local Variables
4453  TYPE(field_type), POINTER :: GEOMETRIC_FIELD
4454  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4455  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
4456  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4457  TYPE(varying_string) :: LOCAL_ERROR
4458 
4459  REAL(DP) :: ALPHA
4460 
4461 
4462  enters("Darcy_PreSolveStorePreviousData",err,error,*999)
4463 
4464  NULLIFY(solver_equations)
4465  NULLIFY(solver_mapping)
4466  NULLIFY(equations_set)
4467 
4468  IF(ASSOCIATED(control_loop)) THEN
4469  IF(ASSOCIATED(solver)) THEN
4470  IF(solver%GLOBAL_NUMBER==solver_number_darcy) THEN
4471  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4472  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4473  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4474  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4475  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
4476  END IF
4477  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4478  CASE(problem_standard_darcy_subtype)
4479  ! do nothing
4480  CASE(problem_quasistatic_darcy_subtype)
4481  ! do nothing
4482  CASE(problem_transient_darcy_subtype)
4483  ! do nothing
4484  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype,problem_standard_elasticity_darcy_subtype, &
4485  & problem_pgm_elasticity_darcy_subtype,problem_pgm_transient_darcy_subtype, &
4486  & problem_quasistatic_elasticity_transient_darcy_subtype,problem_quasistatic_elast_trans_darcy_mat_solve_subtype)
4487  solver_equations=>solver%SOLVER_EQUATIONS
4488  IF(ASSOCIATED(solver_equations)) THEN
4489  solver_mapping=>solver_equations%SOLVER_MAPPING
4490  IF(ASSOCIATED(solver_mapping)) THEN
4491  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
4492  IF(ASSOCIATED(equations_set)) THEN
4493  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
4494  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
4495  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
4496  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
4497  & err,error,*999)
4498  END IF
4499  SELECT CASE(equations_set%SPECIFICATION(3))
4500  CASE(equations_set_standard_darcy_subtype)
4501  ! do nothing
4502  CASE(equations_set_quasistatic_darcy_subtype)
4503  ! do nothing
4504  CASE(equations_set_transient_darcy_subtype)
4505  ! do nothing
4506  CASE(equations_set_ale_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
4507  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
4508  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
4509  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
4510  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
4511  IF(ASSOCIATED(geometric_field)) THEN
4512  !--- Store the GEOMETRY field values of the previous time step
4513  alpha = 1.0_dp
4514  CALL field_parameter_sets_copy(geometric_field,field_u_variable_type, &
4515  & field_values_set_type,field_previous_values_set_type,alpha,err,error,*999)
4516  ELSE
4517  CALL flagerror("Dependent field and / or geometric field is / are not associated.",err,error,*999)
4518  ENDIF
4519  CASE DEFAULT
4520  local_error="Equations set subtype " &
4521  & //trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
4522  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4523  CALL flagerror(local_error,err,error,*999)
4524  END SELECT
4525  ELSE
4526  CALL flagerror("Equations set is not associated.",err,error,*999)
4527  ENDIF
4528  ELSE
4529  CALL flagerror("Solver mapping is not associated.",err,error,*999)
4530  ENDIF
4531  ELSE
4532  CALL flagerror("Solver equations is not associated.",err,error,*999)
4533  ENDIF
4534  CASE DEFAULT
4535  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4536  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4537  CALL flagerror(local_error,err,error,*999)
4538  END SELECT
4539  ELSE
4540  CALL flagerror("Problem is not associated.",err,error,*999)
4541  ENDIF
4542  ELSE
4543  ! do nothing
4544  ENDIF
4545  ELSE
4546  CALL flagerror("Solver is not associated.",err,error,*999)
4547  ENDIF
4548  ELSE
4549  CALL flagerror("Control loop is not associated.",err,error,*999)
4550  ENDIF
4551 
4552  exits("Darcy_PreSolveStorePreviousData")
4553  RETURN
4554 999 errorsexits("Darcy_PreSolveStorePreviousData",err,error)
4555  RETURN 1
4556 
4557  END SUBROUTINE darcy_presolvestorepreviousdata
4558 
4559  !
4560  !================================================================================================================================
4561  !
4562 
4564  SUBROUTINE darcy_equation_pre_solve_ale_update_mesh(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4565 
4566  !Argument variables
4567  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4568  TYPE(solver_type), POINTER :: SOLVER
4569  INTEGER(INTG), INTENT(OUT) :: ERR
4570  TYPE(varying_string), INTENT(OUT) :: ERROR
4571  !Local Variables
4572  TYPE(field_type), POINTER :: GEOMETRIC_FIELD
4573  TYPE(solver_type), POINTER :: SOLVER_ALE_DARCY
4574  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4575  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
4576  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4577  TYPE(control_loop_type), POINTER :: CONTROL_TIME_LOOP
4578  TYPE(varying_string) :: LOCAL_ERROR
4579 
4580  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT,ALPHA
4581  REAL(DP), POINTER :: MESH_DISPLACEMENT_VALUES(:)
4582 
4583  INTEGER(INTG) :: dof_number,NUMBER_OF_DOFS,NDOFS_TO_PRINT,loop_idx
4584  INTEGER(INTG) :: PROBLEM_SUBTYPE
4585 
4586  enters("DARCY_EQUATION_PRE_SOLVE_ALE_UPDATE_MESH",err,error,*999)
4587 
4588  NULLIFY(solver_ale_darcy)
4589  NULLIFY(solver_equations)
4590  NULLIFY(solver_mapping)
4591  NULLIFY(equations_set)
4592 
4593  IF(ASSOCIATED(control_loop)) THEN
4594  control_time_loop=>control_loop
4595  DO loop_idx=1,control_loop%CONTROL_LOOP_LEVEL
4596  IF(control_time_loop%LOOP_TYPE==problem_control_time_loop_type) THEN
4597  CALL control_loop_current_times_get(control_time_loop,current_time,time_increment,err,error,*999)
4598  EXIT
4599  ENDIF
4600  IF (ASSOCIATED(control_loop%PARENT_LOOP)) THEN
4601  control_time_loop=>control_time_loop%PARENT_LOOP
4602  ELSE
4603  CALL flagerror("Could not find a time control loop.",err,error,*999)
4604  ENDIF
4605  ENDDO
4606 
4607  IF(ASSOCIATED(solver)) THEN
4608  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4609  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4610  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4611  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4612  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
4613  END IF
4614  problem_subtype=control_loop%PROBLEM%SPECIFICATION(3)
4615  SELECT CASE(problem_subtype)
4616  CASE(problem_standard_darcy_subtype)
4617  ! do nothing
4618  CASE(problem_quasistatic_darcy_subtype)
4619  ! do nothing
4620  CASE(problem_transient_darcy_subtype)
4621  ! do nothing
4622  CASE(problem_ale_darcy_subtype,problem_standard_elasticity_darcy_subtype,problem_pgm_darcy_subtype, &
4623  & problem_pgm_elasticity_darcy_subtype,problem_pgm_transient_darcy_subtype, &
4624  & problem_quasistatic_elasticity_transient_darcy_subtype,problem_quasistatic_elast_trans_darcy_mat_solve_subtype)
4625  IF(solver%GLOBAL_NUMBER==solver_number_darcy) THEN
4626  solver_equations=>solver%SOLVER_EQUATIONS
4627  IF(ASSOCIATED(solver_equations)) THEN
4628  solver_mapping=>solver_equations%SOLVER_MAPPING
4629  IF(ASSOCIATED(solver_mapping)) THEN
4630  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
4631  IF(ASSOCIATED(equations_set)) THEN
4632  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
4633  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
4634  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
4635  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
4636  & err,error,*999)
4637  END IF
4638  SELECT CASE(equations_set%SPECIFICATION(3))
4639  CASE(equations_set_standard_darcy_subtype)
4640  ! do nothing
4641  CASE(equations_set_quasistatic_darcy_subtype)
4642  ! do nothing
4643  CASE(equations_set_ale_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
4644  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
4645  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
4646  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
4647  IF(solver%OUTPUT_TYPE>=solver_progress_output) THEN
4648  CALL write_string(general_output_type,"Darcy update mesh ... ",err,error,*999)
4649  ENDIF
4650  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
4651  IF(ASSOCIATED(geometric_field)) THEN
4652  !--- First, get pointer to mesh displacement values
4653  NULLIFY(mesh_displacement_values)
4654  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type, &
4655  & field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
4656  IF(diagnostics1) THEN
4657  ndofs_to_print = SIZE(mesh_displacement_values,1)
4658  CALL write_string_vector(diagnostic_output_type,1,1,ndofs_to_print,ndofs_to_print,ndofs_to_print,&
4659  & mesh_displacement_values,'(" MESH_DISPLACEMENT_VALUES = ",3(X,E13.6))','3(3(X,E13.6))', &
4660  & err,error,*999)
4661  ENDIF
4662 
4663  number_of_dofs = geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR%NUMBER_OF_DOFS
4664 
4665  IF(problem_subtype==problem_quasistatic_elasticity_transient_darcy_subtype &
4666  & .OR. problem_subtype==problem_quasistatic_elast_trans_darcy_mat_solve_subtype &
4667  & .OR. problem_subtype==problem_standard_elasticity_darcy_subtype) THEN
4668  !--- Don't update geometric field here, this is done in
4669  ! darcy_equation_pre_solve_get_solid_displacement for these problems, but
4670  ! needs to be made consistent between the different problem types
4671  ELSE
4672  !--- Second, update geometric field
4673  DO dof_number=1,number_of_dofs
4674  CALL field_parameter_set_add_local_dof(geometric_field, &
4675  & field_u_variable_type,field_values_set_type,dof_number, &
4676  & mesh_displacement_values(dof_number), &
4677  & err,error,*999)
4678  END DO
4679  CALL field_parameter_set_update_start(geometric_field, &
4680  & field_u_variable_type, field_values_set_type,err,error,*999)
4681  CALL field_parameter_set_update_finish(geometric_field, &
4682  & field_u_variable_type, field_values_set_type,err,error,*999)
4683  ENDIF
4684 
4685  !--- Third, use displacement values to calculate velocity values
4686  alpha=1.0_dp/time_increment
4687  CALL field_parameter_sets_copy(geometric_field,field_u_variable_type, &
4688  & field_mesh_displacement_set_type,field_mesh_velocity_set_type,alpha,err,error,*999)
4689  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type, &
4690  & field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
4691  ELSE
4692  CALL flagerror("Geometric field is not associated.",err,error,*999)
4693  ENDIF
4694  CASE DEFAULT
4695  local_error="Equations set subtype " &
4696  & //trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
4697  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4698  CALL flagerror(local_error,err,error,*999)
4699  END SELECT
4700  ELSE
4701  CALL flagerror("Equations set is not associated.",err,error,*999)
4702  ENDIF
4703  ELSE
4704  CALL flagerror("Solver mapping is not associated.",err,error,*999)
4705  ENDIF
4706  ELSE
4707  CALL flagerror("Solver equations is not associated.",err,error,*999)
4708  ENDIF
4709  ELSE
4710  ! do nothing
4711  ENDIF
4712  CASE DEFAULT
4713  local_error="Problem subtype "//trim(number_to_vstring(problem_subtype,"*",err,error))// &
4714  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4715  CALL flagerror(local_error,err,error,*999)
4716  END SELECT
4717  ELSE
4718  CALL flagerror("Problem is not associated.",err,error,*999)
4719  ENDIF
4720  ELSE
4721  CALL flagerror("Solver is not associated.",err,error,*999)
4722  ENDIF
4723  ELSE
4724  CALL flagerror("Control loop is not associated.",err,error,*999)
4725  ENDIF
4726 
4727  exits("DARCY_EQUATION_PRE_SOLVE_ALE_UPDATE_MESH")
4728  RETURN
4729 999 errorsexits("DARCY_EQUATION_PRE_SOLVE_ALE_UPDATE_MESH",err,error)
4730  RETURN 1
4731  END SUBROUTINE darcy_equation_pre_solve_ale_update_mesh
4732 
4733  !
4734  !================================================================================================================================
4735  !
4736 
4738  SUBROUTINE darcy_presolveupdateboundaryconditions(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4739 
4740  !Argument variables
4741  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4742  TYPE(solver_type), POINTER :: SOLVER
4743  INTEGER(INTG), INTENT(OUT) :: ERR
4744  TYPE(varying_string), INTENT(OUT) :: ERROR
4745  !Local Variables
4746  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4747  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
4748  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4749  TYPE(equations_type), POINTER :: EQUATIONS
4750  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
4751  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
4752  TYPE(varying_string) :: LOCAL_ERROR
4753  TYPE(field_type), POINTER :: DEPENDENT_FIELD, GEOMETRIC_FIELD
4754  TYPE(boundary_conditions_variable_type), POINTER :: BOUNDARY_CONDITIONS_VARIABLE
4755  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
4756  TYPE(control_loop_type), POINTER :: CONTROL_TIME_LOOP
4757 
4758  REAL(DP), POINTER :: MESH_VELOCITY_VALUES(:)
4759  REAL(DP), POINTER :: INITIAL_VALUES(:)
4760  REAL(DP), POINTER :: DUMMY_VALUES1(:)
4761  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
4762  REAL(DP) :: PRESSURE
4763 
4764  INTEGER(INTG) :: FIELD_VAR_TYPE
4765  INTEGER(INTG) :: BOUNDARY_CONDITION_CHECK_VARIABLE
4766  INTEGER(INTG) :: dof_number,NUMBER_OF_DOFS,loop_idx
4767  INTEGER(INTG) :: NDOFS_TO_PRINT
4768 
4769  enters("Darcy_PreSolveUpdateBoundaryConditions",err,error,*999)
4770 
4771  IF(ASSOCIATED(control_loop)) THEN
4772  control_time_loop=>control_loop
4773  DO loop_idx=1,control_loop%CONTROL_LOOP_LEVEL
4774  IF(control_time_loop%LOOP_TYPE==problem_control_time_loop_type) THEN
4775  CALL control_loop_current_times_get(control_time_loop,current_time,time_increment,err,error,*999)
4776  EXIT
4777  ENDIF
4778  IF (ASSOCIATED(control_loop%PARENT_LOOP)) THEN
4779  control_time_loop=>control_time_loop%PARENT_LOOP
4780  ELSE
4781  CALL flagerror("Could not find a time control loop.",err,error,*999)
4782  ENDIF
4783  ENDDO
4784 
4785  IF(ASSOCIATED(solver)) THEN
4786  IF(solver%GLOBAL_NUMBER==solver_number_darcy) THEN
4787  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4788  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4789  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4790  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4791  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
4792  END IF
4793  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4794  CASE(problem_standard_darcy_subtype)
4795  ! do nothing
4796  CASE(problem_quasistatic_darcy_subtype)
4797  ! do nothing
4798  CASE(problem_transient_darcy_subtype)
4799  ! do nothing
4800  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype,problem_standard_elasticity_darcy_subtype, &
4801  & problem_pgm_elasticity_darcy_subtype,problem_pgm_transient_darcy_subtype, &
4802  & problem_quasistatic_elasticity_transient_darcy_subtype,problem_quasistatic_elast_trans_darcy_mat_solve_subtype)
4803  solver_equations=>solver%SOLVER_EQUATIONS
4804  IF(ASSOCIATED(solver_equations)) THEN
4805  solver_mapping=>solver_equations%SOLVER_MAPPING
4806  IF(ASSOCIATED(solver_mapping)) THEN
4807  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
4808  IF(ASSOCIATED(equations)) THEN
4809  equations_set=>equations%EQUATIONS_SET
4810  IF(ASSOCIATED(equations_set)) THEN
4811  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
4812  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
4813  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
4814  CALL flagerror("Equations set specification must have three entries for a Darcy type equations set.", &
4815  & err,error,*999)
4816  END IF
4817  SELECT CASE(equations_set%SPECIFICATION(3))
4818  CASE(equations_set_standard_darcy_subtype)
4819  ! do nothing
4820  CASE(equations_set_quasistatic_darcy_subtype)
4821  ! do nothing
4822  CASE(equations_set_ale_darcy_subtype,equations_set_incompressible_finite_elasticity_darcy_subtype, &
4823  & equations_set_transient_ale_darcy_subtype,equations_set_elasticity_darcy_inria_model_subtype, &
4824  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
4825  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
4826  IF(solver%OUTPUT_TYPE>=solver_progress_output) THEN
4827  CALL write_string(general_output_type,"Darcy update boundary conditions ... ",err,error,*999)
4828  ENDIF
4829  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
4830  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
4831  IF(ASSOCIATED(dependent_field).AND.ASSOCIATED(geometric_field)) THEN
4832  boundary_conditions=>solver_equations%BOUNDARY_CONDITIONS
4833  IF(ASSOCIATED(boundary_conditions)) THEN
4834  equations_mapping=>equations_set%EQUATIONS%EQUATIONS_MAPPING
4835  IF(ASSOCIATED(equations_mapping)) THEN
4836 
4837  SELECT CASE(equations_set%SPECIFICATION(3))
4838  CASE(equations_set_ale_darcy_subtype, &
4839  & equations_set_incompressible_finite_elasticity_darcy_subtype)
4840  field_variable=>equations_mapping%LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
4841  ! '1' associated with linear matrix
4842  CASE(equations_set_transient_ale_darcy_subtype, &
4843  & equations_set_elasticity_darcy_inria_model_subtype, &
4844  & equations_set_incompressible_elasticity_driven_darcy_subtype, &
4845  & equations_set_incompressible_elast_multi_comp_darcy_subtype)
4846  field_variable=>equations_mapping%DYNAMIC_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
4847  END SELECT
4848 
4849  IF(ASSOCIATED(field_variable)) THEN
4850  field_var_type=field_variable%VARIABLE_TYPE
4851 
4852  CALL boundary_conditions_variable_get(boundary_conditions,field_variable, &
4853  & boundary_conditions_variable,err,error,*999)
4854  IF(ASSOCIATED(boundary_conditions_variable)) THEN
4855  NULLIFY(mesh_velocity_values)
4856  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type, &
4857  & field_mesh_velocity_set_type,mesh_velocity_values,err,error,*999)
4858  NULLIFY(initial_values)
4859  CALL field_parameter_set_data_get(dependent_field,field_var_type, &
4860  & field_initial_values_set_type,initial_values,err,error,*999)
4861  IF(diagnostics1) THEN
4862  NULLIFY( dummy_values1 )
4863  CALL field_parameter_set_data_get(dependent_field,field_var_type, &
4864  & field_values_set_type,dummy_values1,err,error,*999)
4865  ndofs_to_print = SIZE(dummy_values1,1)
4866  CALL write_string_vector(diagnostic_output_type,1,1,ndofs_to_print,ndofs_to_print, &
4867  & ndofs_to_print,dummy_values1, &
4868  & '(" DEPENDENT_FIELD,FIELD_VAR_TYPE,FIELD_VALUES_SET_TYPE (before) = ",4(X,E13.6))', &
4869  & '4(4(X,E13.6))',err,error,*999)
4870  ENDIF
4871  number_of_dofs = dependent_field%VARIABLE_TYPE_MAP(field_var_type)%PTR%NUMBER_OF_DOFS
4872  DO dof_number=1,number_of_dofs
4873  boundary_condition_check_variable=boundary_conditions_variable% &
4874  & condition_types(dof_number)
4875  IF(boundary_condition_check_variable==boundary_condition_moved_wall) THEN
4876  !--- Reset boundary condition to the initial normal-velocity boundary condition
4877  CALL field_parameter_set_update_local_dof(dependent_field, &
4878  & field_var_type,field_values_set_type,dof_number, &
4879  & initial_values(dof_number),err,error,*999)
4880  !--- Add the velocity of the moving boundary on top of the initial boundary condition
4881  !! === If we solve in terms of Darcy flow vector, then do not add mesh velocity === !!
4882  !! === The BC is kept to the initial BC, for instance: null-flux === !!
4883 ! CALL FIELD_PARAMETER_SET_ADD_LOCAL_DOF(DEPENDENT_FIELD, &
4884 ! & FIELD_VAR_TYPE,FIELD_VALUES_SET_TYPE,dof_number, &
4885 ! & MESH_VELOCITY_VALUES(dof_number),ERR,ERROR,*999)
4886 ! ! dependent field ( V_u, V_v, V_w, P_p )
4887 ! ! MESH_VELOCITY_VALUES ( V_u, V_v, V_w )
4888 
4889  ELSE IF( boundary_condition_check_variable==boundary_condition_fixed .AND. &
4890  & equations_set%SPECIFICATION(3)==equations_set_elasticity_darcy_inria_model_subtype) THEN
4891  !\ToDo: Check component number; this way we can also apply it to velocity
4892  !--- Set the time-dependent pressure BC
4893  pressure = initial_values(dof_number) * (1.0_dp - exp(- current_time**2.0_dp / 0.25_dp))
4894 
4895  CALL field_parameter_set_update_local_dof(dependent_field, &
4896  & field_var_type,field_values_set_type,dof_number, &
4897  & pressure,err,error,*999)
4898  ELSE
4899  ! do nothing
4900  END IF
4901  END DO
4902  CALL field_parameter_set_update_start(dependent_field, &
4903  & field_var_type, field_values_set_type,err,error,*999)
4904  CALL field_parameter_set_update_finish(dependent_field, &
4905  & field_var_type, field_values_set_type,err,error,*999)
4906  IF(diagnostics1) THEN
4907  ndofs_to_print = SIZE(mesh_velocity_values,1)
4908  CALL write_string_vector(diagnostic_output_type,1,1,ndofs_to_print,ndofs_to_print, &
4909  & ndofs_to_print,mesh_velocity_values, &
4910  & '(" MESH_VELOCITY_VALUES = ",4(X,E13.6))','4(4(X,E13.6))',err,error,*999)
4911  CALL write_string(diagnostic_output_type," ",err,error,*999)
4912  !
4913  ndofs_to_print = SIZE(initial_values,1)
4914  CALL write_string_vector(diagnostic_output_type,1,1,ndofs_to_print,ndofs_to_print, &
4915  & ndofs_to_print,initial_values, &
4916  & '(" INITIAL_VALUES = ",4(X,E13.6))', &
4917  & '4(4(X,E13.6))',err,error,*999)
4918  !
4919  NULLIFY( dummy_values1 )
4920  CALL field_parameter_set_data_get(dependent_field,field_var_type, &
4921  & field_values_set_type,dummy_values1,err,error,*999)
4922  ndofs_to_print = SIZE(dummy_values1,1)
4923  CALL write_string_vector(diagnostic_output_type,1,1,ndofs_to_print,ndofs_to_print, &
4924  & ndofs_to_print,dummy_values1, &
4925  & '(" DEPENDENT_FIELD,FIELD_VAR_TYPE,FIELD_VALUES_SET_TYPE (after) = ",4(X,E13.6))', &
4926  & '4(4(X,E13.6))',err,error,*999)
4927  CALL field_parameter_set_data_restore(dependent_field,field_var_type, &
4928  & field_values_set_type,dummy_values1,err,error,*999)
4929  ENDIF
4930  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type, &
4931  & field_mesh_velocity_set_type,mesh_velocity_values,err,error,*999)
4932  CALL field_parameter_set_data_restore(dependent_field,field_var_type, &
4933  & field_initial_values_set_type,initial_values,err,error,*999)
4934  ELSE
4935  CALL flagerror("Boundary condition variable is not associated.",err,error,*999)
4936  END IF
4937 
4938  CALL field_parameter_set_update_start(equations_set%DEPENDENT%DEPENDENT_FIELD,field_var_type, &
4939  & field_values_set_type,err,error,*999)
4940  CALL field_parameter_set_update_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,field_var_type, &
4941  & field_values_set_type,err,error,*999)
4942 
4943  ELSE
4944  CALL flagerror("FIELD_VAR_TYPE is not associated.",err,error,*999)
4945  ENDIF
4946  ELSE
4947  CALL flagerror("EQUATIONS_MAPPING is not associated.",err,error,*999)
4948  ENDIF
4949  ELSE
4950  CALL flagerror("Boundary conditions are not associated.",err,error,*999)
4951  END IF
4952  ELSE
4953  CALL flagerror("Dependent field and/or geometric field is/are not associated.",err,error,*999)
4954  END IF
4955  CASE DEFAULT
4956  local_error="Equations set subtype " &
4957  & //trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
4958  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4959  CALL flagerror(local_error,err,error,*999)
4960  END SELECT
4961  ELSE
4962  CALL flagerror("Equations set is not associated.",err,error,*999)
4963  END IF
4964  ELSE
4965  CALL flagerror("Equations are not associated.",err,error,*999)
4966  END IF
4967  ELSE
4968  CALL flagerror("Solver mapping is not associated.",err,error,*999)
4969  ENDIF
4970  ELSE
4971  CALL flagerror("Solver equations are not associated.",err,error,*999)
4972  END IF
4973  CASE DEFAULT
4974  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4975  & " is not valid for a Darcy equation fluid type of a fluid mechanics problem class."
4976  CALL flagerror(local_error,err,error,*999)
4977  END SELECT
4978  ELSE
4979  CALL flagerror("Problem is not associated.",err,error,*999)
4980  ENDIF
4981  ELSE
4982  ! do nothing
4983  ENDIF
4984  ELSE
4985  CALL flagerror("Solver is not associated.",err,error,*999)
4986  ENDIF
4987  ELSE
4988  CALL flagerror("Control loop is not associated.",err,error,*999)
4989  ENDIF
4990 
4991  exits("Darcy_PreSolveUpdateBoundaryConditions")
4992  RETURN
4993 999 errors("Darcy_PreSolveUpdateBoundaryConditions",err,error)
4994  exits("Darcy_PreSolveUpdateBoundaryConditions")
4995  RETURN 1
4996 
4997  END SUBROUTINE darcy_presolveupdateboundaryconditions
4998 
4999  !
5000  !================================================================================================================================
5001  !
5002 
5004  SUBROUTINE darcy_presolveupdatematrixproperties(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
5005 
5006  !Argument variables
5007  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
5008  TYPE(solver_type), POINTER :: SOLVER
5009  INTEGER(INTG), INTENT(OUT) :: ERR
5010  TYPE(varying_string), INTENT(OUT) :: ERROR
5011 
5012  !Local Variables
5013  TYPE(solver_type), POINTER :: SOLVER_MAT_PROPERTIES, SOLVER_ALE_DARCY
5014  TYPE(field_type), POINTER :: DEPENDENT_FIELD_MAT_PROPERTIES, MATERIALS_FIELD_ALE_DARCY
5015  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS_MAT_PROPERTIES, SOLVER_EQUATIONS_ALE_DARCY
5016  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING_MAT_PROPERTIES, SOLVER_MAPPING_ALE_DARCY
5017  TYPE(equations_set_type), POINTER :: EQUATIONS_SET_MAT_PROPERTIES, EQUATIONS_SET_ALE_DARCY
5018  TYPE(varying_string) :: LOCAL_ERROR
5019 
5020  REAL(DP), POINTER :: DUMMY_VALUES2(:)
5021 
5022  INTEGER(INTG) :: NUMBER_OF_COMPONENTS_DEPENDENT_FIELD_MAT_PROPERTIES,NUMBER_OF_COMPONENTS_MATERIALS_FIELD_ALE_DARCY
5023  INTEGER(INTG) :: NDOFS_TO_PRINT, I
5024 
5025 
5026  enters("Darcy_PreSolveUpdateMatrixProperties",err,error,*999)
5027 
5028  IF(ASSOCIATED(control_loop)) THEN
5029 
5030  NULLIFY(solver_mat_properties)
5031  NULLIFY(solver_ale_darcy)
5032 
5033  IF(ASSOCIATED(solver)) THEN
5034  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
5035  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
5036  CALL flagerror("Problem specification is not allocated.",err,error,*999)
5037  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
5038  CALL flagerror("Problem specification must have three entries for a Darcy equation problem.",err,error,*999)
5039  END IF
5040  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
5041  CASE(problem_standard_darcy_subtype)
5042  ! do nothing
5043  CASE(problem_quasistatic_darcy_subtype)
5044  ! do nothing
5045  CASE(problem_transient_darcy_subtype)
5046  ! do nothing
5047  CASE(problem_ale_darcy_subtype,problem_pgm_darcy_subtype,problem_standard_elasticity_darcy_subtype, &