OpenCMISS-Iron Internal API Documentation
Poisson_equations_routines.f90
Go to the documentation of this file.
1 
43 
46 
47  USE base_routines
48  USE basis_routines
50  USE constants
54  USE domain_mappings
59  USE field_routines
60  USE input_output
62  USE kinds
63  USE maths
64  USE matrix_vector
65  USE node_routines
67  USE strings
68  USE solver_routines
69  USE timer
70  USE types
71 ! temporary input for vector-source
73 
74 #include "macros.h"
75 
76  IMPLICIT NONE
77 
78  PRIVATE
79 
80  !Module parameters
81 
82  !Module types
83 
84  !Module variables
85 
86  !Interfaces
87 
89 
91 
93 
95 
97 
99 
101 
103 
105 
106 CONTAINS
107 
108  !
109  !================================================================================================================================
110  !
111 
113  SUBROUTINE poisson_boundaryconditionsanalyticcalculate(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*)
115  !Argument variables
116  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
117  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
118  INTEGER(INTG), INTENT(OUT) :: ERR
119  TYPE(varying_string), INTENT(OUT) :: ERROR
120  !Local Variables
121  INTEGER(INTG) :: component_idx,deriv_idx,dim_idx,local_ny,node_idx,NUMBER_OF_DIMENSIONS,variable_idx,variable_type,global_ny
122  REAL(DP) :: VALUE,X(3)
123  REAL(DP), POINTER :: GEOMETRIC_PARAMETERS(:)
124  TYPE(domain_type), POINTER :: DOMAIN
125  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
126  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
127  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
128  TYPE(varying_string) :: LOCAL_ERROR
129 
130  enters("Poisson_BoundaryConditionsAnalyticCalculate",err,error,*999)
131 
132  IF(ASSOCIATED(equations_set)) THEN
133  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
134  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
135  IF(ASSOCIATED(dependent_field)) THEN
136  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
137  IF(ASSOCIATED(geometric_field)) THEN
138  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
139  NULLIFY(geometric_variable)
140  CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
141  NULLIFY(geometric_parameters)
142  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,geometric_parameters, &
143  & err,error,*999)
144 
145  IF(ASSOCIATED(boundary_conditions)) THEN
146 
147  DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES !U and deludeln
148  variable_type=dependent_field%VARIABLES(variable_idx)%VARIABLE_TYPE
149  field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
150  IF(ASSOCIATED(field_variable)) THEN
151  CALL field_parameter_set_create(dependent_field,variable_type,field_analytic_values_set_type,err,error,*999)
152  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS !u,v,w
153  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
154  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
155  IF(ASSOCIATED(domain)) THEN
156  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
157  domain_nodes=>domain%TOPOLOGY%NODES
158  IF(ASSOCIATED(domain_nodes)) THEN
159 
160  !Loop over the local nodes excluding the ghosts.
161  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
162  !!TODO \todo We should interpolate the geometric field here and the node position.
163  DO dim_idx=1,number_of_dimensions
164  local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
165  & nodes(node_idx)%DERIVATIVES(1)%VERSIONS(1)
166  x(dim_idx)=geometric_parameters(local_ny)
167  ENDDO !dim_idx
168  !Loop over the derivatives
169  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
170  SELECT CASE(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE)
172  !u=ln(4/(x+y+1^2))
173  SELECT CASE(variable_type)
174  CASE(field_u_variable_type)
175  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
176  CASE(no_global_deriv)
177  VALUE=log(4.0_dp/((x(1)+x(2)+1.0_dp)**2))
178  CASE(global_deriv_s1)
179  CALL flagerror("Not implemented.",err,error,*999)
180  CASE(global_deriv_s2)
181  CALL flagerror("Not implemented.",err,error,*999)
182  CASE(global_deriv_s1_s2)
183  CALL flagerror("Not implemented.",err,error,*999)
184  CASE DEFAULT
185  local_error="The global derivative index of "//trim(number_to_vstring( &
186  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
187  & err,error))//" is invalid."
188  CALL flagerror(local_error,err,error,*999)
189  END SELECT
190  CASE(field_deludeln_variable_type)
191  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
192  CASE(no_global_deriv)
193  !This is believed to be incorrect, should be: VALUE=-2.0_DP/(X(1)+X(2)+1.0_DP)
194  VALUE=-2.0_dp*(x(1)+x(2))/(x(1)+x(2)+1.0_dp)
195  CASE(global_deriv_s1)
196  CALL flagerror("Not implemented.",err,error,*999)
197  CASE(global_deriv_s2)
198  CALL flagerror("Not implemented.",err,error,*999)
199  CASE(global_deriv_s1_s2)
200  CALL flagerror("Not implemented.",err,error,*999)
201  CASE DEFAULT
202  local_error="The global derivative index of "//trim(number_to_vstring( &
203  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
204  & err,error))//" is invalid."
205  CALL flagerror(local_error,err,error,*999)
206  END SELECT
207  END SELECT
209  CALL flagerror("The analytic function type is not implemented yet.",err,error,*999)
211  CALL flagerror("The analytic function type is not implemented yet.",err,error,*999)
213  !u=ln(6/(x+y+z+1^2))
214  SELECT CASE(variable_type)
215  CASE(field_u_variable_type)
216  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
217  CASE(no_global_deriv)
218  VALUE=log(6.0_dp/((x(1)+x(2)+x(3)+1.0_dp)**2))
219  CASE(global_deriv_s1)
220  CALL flagerror("Not implemented.",err,error,*999)
221  CASE(global_deriv_s2)
222  CALL flagerror("Not implemented.",err,error,*999)
223  CASE(global_deriv_s1_s2)
224  CALL flagerror("Not implemented.",err,error,*999)
225  CASE DEFAULT
226  local_error="The global derivative index of "//trim(number_to_vstring( &
227  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
228  & err,error))//" is invalid."
229  CALL flagerror(local_error,err,error,*999)
230  END SELECT
231  CASE(field_deludeln_variable_type)
232  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
233  CASE(no_global_deriv)
234  VALUE=-3.0_dp/(x(1)+x(2)+x(3)+1.0_dp)
235  CASE(global_deriv_s1)
236  CALL flagerror("Not implemented.",err,error,*999)
237  CASE(global_deriv_s2)
238  CALL flagerror("Not implemented.",err,error,*999)
239  CASE(global_deriv_s1_s2)
240  CALL flagerror("Not implemented.",err,error,*999)
241  CASE DEFAULT
242  local_error="The global derivative index of "//trim(number_to_vstring( &
243  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
244  & err,error))//" is invalid."
245  CALL flagerror(local_error,err,error,*999)
246  END SELECT
247  CASE DEFAULT
248  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
249  & " is invalid."
250  CALL flagerror(local_error,err,error,*999)
251  END SELECT
253  CALL flagerror("The analytic function type is not implemented yet.",err,error,*999)
255  CALL flagerror("The analytic function type is not implemented yet.",err,error,*999)
257  !\todo: This test case has been set up for the Pressure Poisson equation - needs to be formulated in a more general way
258  !u=ln(4/(x+y+1^2))
259  SELECT CASE(variable_type)
260  CASE(field_u_variable_type)
261  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
262  CASE(no_global_deriv)
263  VALUE=sin(2.0_dp*pi*x(1)/10.0_dp)*sin(2.0_dp*pi*x(2)/10.0_dp)*sin(2.0_dp*pi*x(3)/10.0_dp)
264  ! VALUE=2*X(1)*X(1)+2*X(2)
265  CASE(global_deriv_s1)
266  CALL flagerror("Not implemented.",err,error,*999)
267  CASE(global_deriv_s2)
268  CALL flagerror("Not implemented.",err,error,*999)
269  CASE(global_deriv_s1_s2)
270  CALL flagerror("Not implemented.",err,error,*999)
271  CASE DEFAULT
272  local_error="The global derivative index of "//trim(number_to_vstring( &
273  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
274  & err,error))//" is invalid."
275  CALL flagerror(local_error,err,error,*999)
276  END SELECT
277  CASE(field_deludeln_variable_type)
278  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
279  CASE(no_global_deriv)
280  ! do nothing, tbd
281  CASE(global_deriv_s1)
282  CALL flagerror("Not implemented.",err,error,*999)
283  CASE(global_deriv_s2)
284  CALL flagerror("Not implemented.",err,error,*999)
285  CASE(global_deriv_s1_s2)
286  CALL flagerror("Not implemented.",err,error,*999)
287  CASE DEFAULT
288  local_error="The global derivative index of "//trim(number_to_vstring( &
289  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
290  & err,error))//" is invalid."
291  CALL flagerror(local_error,err,error,*999)
292  END SELECT
293  CASE DEFAULT
294  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
295  & " is invalid."
296  CALL flagerror(local_error,err,error,*999)
297  END SELECT
298  CASE DEFAULT
299  local_error="The analytic function type of "// &
300  & trim(number_to_vstring(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
301  & " is invalid."
302  CALL flagerror(local_error,err,error,*999)
303  END SELECT
304  !Default to version 1 of each node derivative
305  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
306  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
307  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
308  & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
309 
310  global_ny=field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
311 
312  IF(variable_type==field_u_variable_type.AND.domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN !For Dirichlet
313  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field,variable_type, &
314  & local_ny,boundary_condition_fixed,VALUE,err,error,*999)
315  ENDIF
316 
317  IF(variable_type==field_deludeln_variable_type.and.node_idx/=1) THEN
318  IF(domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN
319  !If we are a boundary node then set the analytic value on the boundary
320  !CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,DEPENDENT_FIELD,variable_type,local_ny, &
321  ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
322  !Do nothing at present
323  ENDIF
324  ENDIF
325 
326  ENDDO !deriv_idx
327  ENDDO !node_idx
328  ELSE
329  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
330  ENDIF
331  ELSE
332  CALL flagerror("Domain topology is not associated.",err,error,*999)
333  ENDIF
334  ELSE
335  CALL flagerror("Domain is not associated.",err,error,*999)
336  ENDIF
337  ELSE
338  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
339  ENDIF
340  ENDDO !component_idx
341  CALL field_parameter_set_update_start(dependent_field,variable_type,field_analytic_values_set_type, &
342  & err,error,*999)
343  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_analytic_values_set_type, &
344  & err,error,*999)
345  ELSE
346  CALL flagerror("Field variable is not associated.",err,error,*999)
347  ENDIF
348 
349  ENDDO !variable_idx
350 
351  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
352  & geometric_parameters,err,error,*999)
353 
354  ELSE
355  CALL flagerror("Boundary conditions is not associated.",err,error,*999)
356  ENDIF
357  ELSE
358  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
359  ENDIF
360  ELSE
361  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
362  ENDIF
363  ELSE
364  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
365  ENDIF
366  ELSE
367  CALL flagerror("Equations set is not associated.",err,error,*999)
368  ENDIF
369 
370  exits("Poisson_BoundaryConditionsAnalyticCalculate")
371  RETURN
372 999 errorsexits("Poisson_BoundaryConditionsAnalyticCalculate",err,error)
373  RETURN 1
374 
376 
377  !
378  !================================================================================================================================
379  !
380 
382  SUBROUTINE poisson_equation_equations_set_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
384  !Argument variables
385  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
386  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
387  INTEGER(INTG), INTENT(OUT) :: ERR
388  TYPE(varying_string), INTENT(OUT) :: ERROR
389  !Local Variables
390  TYPE(varying_string) :: LOCAL_ERROR
391 
392  enters("POISSON_EQUATION_EQUATIONS_SET_SETUP",err,error,*999)
393 
394  IF(ASSOCIATED(equations_set)) THEN
395  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
396  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
397  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
398  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
399  & err,error,*999)
400  END IF
401  SELECT CASE(equations_set%SPECIFICATION(3))
403  CALL poisson_equationssetlinearsourcesetup(equations_set,equations_set_setup,err,error,*999)
405  CALL poisson_equationssetextracellularbidomainsetup(equations_set,equations_set_setup,err,error,*999)
407  CALL poisson_equationssetlinearsourcesetup(equations_set,equations_set_setup,err,error,*999)
410  CALL poisson_equationssetpressurepoissonsetup(equations_set,equations_set_setup,err,error,*999)
412  CALL poisson_equationssetnonlinearsourcesetup(equations_set,equations_set_setup,err,error,*999)
414  CALL poisson_equationssetnonlinearsourcesetup(equations_set,equations_set_setup,err,error,*999)
415  CASE DEFAULT
416  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
417  & " is not valid for a Poisson equation type of a classical field equation set class."
418  CALL flagerror(local_error,err,error,*999)
419  END SELECT
420  ELSE
421  CALL flagerror("Equations set is not associated.",err,error,*999)
422  ENDIF
423 
424  exits("POISSON_EQUATION_EQUATIONS_SET_SETUP")
425  RETURN
426 999 errorsexits("POISSON_EQUATION_EQUATIONS_SET_SETUP",err,error)
427  RETURN 1
429 
430  !
431  !================================================================================================================================
432  !
433 
435  SUBROUTINE poisson_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)
437  !Argument variables
438  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
439  INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
440  INTEGER(INTG), INTENT(OUT) :: ERR
441  TYPE(varying_string), INTENT(OUT) :: ERROR
442  !Local Variables
443  TYPE(varying_string) :: LOCAL_ERROR
444 
445  enters("Poisson_EquationsSetSolutionMethodSet",err,error,*999)
446 
447  IF(ASSOCIATED(equations_set)) THEN
448  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
449  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
450  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
451  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
452  & err,error,*999)
453  END IF
454  SELECT CASE(equations_set%SPECIFICATION(3))
456  SELECT CASE(solution_method)
458  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
460  CALL flagerror("Not implemented.",err,error,*999)
462  CALL flagerror("Not implemented.",err,error,*999)
464  CALL flagerror("Not implemented.",err,error,*999)
466  CALL flagerror("Not implemented.",err,error,*999)
468  CALL flagerror("Not implemented.",err,error,*999)
469  CASE DEFAULT
470  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
471  CALL flagerror(local_error,err,error,*999)
472  END SELECT
474  SELECT CASE(solution_method)
476  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
478  CALL flagerror("Not implemented.",err,error,*999)
480  CALL flagerror("Not implemented.",err,error,*999)
482  CALL flagerror("Not implemented.",err,error,*999)
484  CALL flagerror("Not implemented.",err,error,*999)
486  CALL flagerror("Not implemented.",err,error,*999)
487  CASE DEFAULT
488  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
489  CALL flagerror(local_error,err,error,*999)
490  END SELECT
494  SELECT CASE(solution_method)
496  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
498  CALL flagerror("Not implemented.",err,error,*999)
500  CALL flagerror("Not implemented.",err,error,*999)
502  CALL flagerror("Not implemented.",err,error,*999)
504  CALL flagerror("Not implemented.",err,error,*999)
506  CALL flagerror("Not implemented.",err,error,*999)
507  CASE DEFAULT
508  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
509  CALL flagerror(local_error,err,error,*999)
510  END SELECT
512  SELECT CASE(solution_method)
514  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
516  CALL flagerror("Not implemented.",err,error,*999)
518  CALL flagerror("Not implemented.",err,error,*999)
520  CALL flagerror("Not implemented.",err,error,*999)
522  CALL flagerror("Not implemented.",err,error,*999)
524  CALL flagerror("Not implemented.",err,error,*999)
525  CASE DEFAULT
526  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
527  CALL flagerror(local_error,err,error,*999)
528  END SELECT
530  SELECT CASE(solution_method)
532  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
534  CALL flagerror("Not implemented.",err,error,*999)
536  CALL flagerror("Not implemented.",err,error,*999)
538  CALL flagerror("Not implemented.",err,error,*999)
540  CALL flagerror("Not implemented.",err,error,*999)
542  CALL flagerror("Not implemented.",err,error,*999)
543  CASE DEFAULT
544  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
545  CALL flagerror(local_error,err,error,*999)
546  END SELECT
547  CASE DEFAULT
548  local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
549  & " is not valid for a Poisson equation type of an classical field equations set class."
550  CALL flagerror(local_error,err,error,*999)
551  END SELECT
552  ELSE
553  CALL flagerror("Equations set is not associated.",err,error,*999)
554  ENDIF
555 
556  exits("Poisson_EquationsSetSolutionMethodSet")
557  RETURN
558 999 errors("Poisson_EquationsSetSolutionMethodSet",err,error)
559  exits("Poisson_EquationsSetSolutionMethodSet")
560  RETURN 1
562 
563  !
564  !================================================================================================================================
565  !
566 
568  SUBROUTINE poisson_equationssetspecificationset(equationsSet,specification,err,error,*)
570  !Argument variables
571  TYPE(equations_set_type), POINTER :: equationsSet
572  INTEGER(INTG), INTENT(IN) :: specification(:)
573  INTEGER(INTG), INTENT(OUT) :: err
574  TYPE(varying_string), INTENT(OUT) :: error
575  !Local Variables
576  TYPE(varying_string) :: localError
577  INTEGER(INTG) :: subtype
578 
579  enters("Poisson_EquationsSetSpecificationSet",err,error,*999)
580 
581  IF(ASSOCIATED(equationsset)) THEN
582  IF(SIZE(specification,1)<3) THEN
583  CALL flagerror("Equations set specification must have at least 3 entries for a Poisson equations set.", &
584  & err,error,*999)
585  END IF
586  subtype=specification(3)
587  SELECT CASE(subtype)
597  !ok
598  CASE DEFAULT
599  localerror="The third equations set specification of "//trim(numbertovstring(subtype,"*",err,error))// &
600  & " is not valid for a Poisson equations set."
601  CALL flagerror(localerror,err,error,*999)
602  END SELECT
603  !Set full specification
604  IF(ALLOCATED(equationsset%specification)) THEN
605  CALL flagerror("Equations set specification is already allocated.",err,error,*999)
606  ELSE
607  ALLOCATE(equationsset%specification(3),stat=err)
608  IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
609  END IF
610  equationsset%specification(1:3)=[equations_set_classical_field_class,equations_set_poisson_equation_type,subtype]
611  ELSE
612  CALL flagerror("Equations set is not associated.",err,error,*999)
613  END IF
614 
615  exits("Poisson_EquationsSetSpecificationSet")
616  RETURN
617 999 errors("Poisson_EquationsSetSpecificationSet",err,error)
618  exits("Poisson_EquationsSetSpecificationSet")
619  RETURN 1
620 
622 
623  !
624  !================================================================================================================================
625  !
626 
628  SUBROUTINE poisson_equationssetpressurepoissonsetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
630  !Argument variables
631  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
632  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
633  INTEGER(INTG), INTENT(OUT) :: ERR
634  TYPE(varying_string), INTENT(OUT) :: ERROR
635  !Local Variables
636  INTEGER(INTG) :: component_idx,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE,NUMBER_OF_DIMENSIONS, &
637  & NUMBER_OF_MATERIALS_COMPONENTS,GEOMETRIC_MESH_COMPONENT,SOURCE_FIELD_NUMBER_OF_COMPONENTS,I, &
638  & SOURCE_FIELD_NUMBER_OF_VARIABLES
639  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
640  TYPE(equations_type), POINTER :: EQUATIONS
641  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
642  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
643  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
644  TYPE(field_type), POINTER :: analytic_field,DEPENDENT_FIELD,geometric_field
645  TYPE(varying_string) :: LOCAL_ERROR
646 
647  enters("Poisson_EquationsSetPressurePoissonSetup",err,error,*999)
648  NULLIFY(equations)
649  NULLIFY(equations_mapping)
650  NULLIFY(equations_matrices)
651  NULLIFY(geometric_decomposition)
652  IF(ASSOCIATED(equations_set)) THEN
653  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
654  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
655  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
656  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
657  & err,error,*999)
658  END IF
659  IF(equations_set%SPECIFICATION(3)==equations_set_linear_pressure_poisson_subtype.OR. &
660  & equations_set%SPECIFICATION(3)==equations_set_nonlinear_pressure_poisson_subtype.OR. &
661  & equations_set%SPECIFICATION(3)==equations_set_ale_pressure_poisson_subtype.OR. &
662  & equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype) THEN
663  SELECT CASE(equations_set_setup%SETUP_TYPE)
665  SELECT CASE(equations_set_setup%ACTION_TYPE)
669 !!TODO: Check valid setup
670  CASE DEFAULT
671  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
672  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
673  & " is invalid for a velocity source Poisson equation."
674  CALL flagerror(local_error,err,error,*999)
675  END SELECT
677  !Do nothing???
679  SELECT CASE(equations_set_setup%ACTION_TYPE)
681  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
682  !Create the auto created dependent field
683  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
684  & dependent_field,err,error,*999)
685  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
686  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
687  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
688  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
689  & err,error,*999)
690  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
691  & geometric_field,err,error,*999)
692  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
693  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,(/field_u_variable_type, &
694  & field_deludeln_variable_type/),err,error,*999)
695  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
696  & field_scalar_dimension_type,err,error,*999)
697  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
698  & field_scalar_dimension_type,err,error,*999)
699  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
700  & field_dp_type,err,error,*999)
701  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
702  & field_dp_type,err,error,*999)
703  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
704  & err,error,*999)
705  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
706  & err,error,*999)
707  !Default to the geometric interpolation setup
708  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
709  & geometric_component_number,err,error,*999)
710  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
711  & geometric_component_number,err,error,*999)
712  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
713  & geometric_component_number,err,error,*999)
714  SELECT CASE(equations_set%SOLUTION_METHOD)
716  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
717  & field_node_based_interpolation,err,error,*999)
718  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
719  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
720  !Default the scaling to the geometric field scaling
721  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
722  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
724  CALL flagerror("Not implemented.",err,error,*999)
726  CALL flagerror("Not implemented.",err,error,*999)
728  CALL flagerror("Not implemented.",err,error,*999)
730  CALL flagerror("Not implemented.",err,error,*999)
732  CALL flagerror("Not implemented.",err,error,*999)
733  CASE DEFAULT
734  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
735  & " is invalid."
736  CALL flagerror(local_error,err,error,*999)
737  END SELECT
738  ELSE
739  !Check the user specified field
740  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
741  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
742  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
743  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type,field_deludeln_variable_type/), &
744  & err,error,*999)
745  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type,err,error,*999)
746  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
747  & err,error,*999)
748  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
749  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
750  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
751  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1,err,error,*999)
752  SELECT CASE(equations_set%SOLUTION_METHOD)
754  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
755  & field_node_based_interpolation,err,error,*999)
756  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
757  & field_node_based_interpolation,err,error,*999)
759  CALL flagerror("Not implemented.",err,error,*999)
761  CALL flagerror("Not implemented.",err,error,*999)
763  CALL flagerror("Not implemented.",err,error,*999)
765  CALL flagerror("Not implemented.",err,error,*999)
767  CALL flagerror("Not implemented.",err,error,*999)
768  CASE DEFAULT
769  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
770  & " is invalid."
771  CALL flagerror(local_error,err,error,*999)
772  END SELECT
773  ENDIF
775  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
776  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
777  !These 3 parameter sets will contain the original velocity field with x,y,z components
778  CALL field_parameter_set_create(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
779  & field_input_vel1_set_type,err,error,*999)
780  CALL field_parameter_set_create(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
781  & field_input_vel2_set_type,err,error,*999)
782  CALL field_parameter_set_create(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
783  & field_input_vel3_set_type,err,error,*999)
784  !This parameter set will contain the inside/outside labelling information for the PPE approach
785  CALL field_parameter_set_create(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
786  & field_input_label_set_type,err,error,*999)
787  ENDIF
788  CASE DEFAULT
789  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
790  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
791  & " is invalid for a velocity source Poisson equation"
792  CALL flagerror(local_error,err,error,*999)
793  END SELECT
795  SELECT CASE(equations_set_setup%ACTION_TYPE)
797  equations_materials=>equations_set%MATERIALS
798  IF(ASSOCIATED(equations_materials)) THEN
799  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
800  !Create the auto created materials field
801  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
802  & materials_field,err,error,*999)
803  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
804  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
805  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
806  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
807  & err,error,*999)
808  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
809  & geometric_field,err,error,*999)
810  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
811  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,(/field_u_variable_type/), &
812  & err,error,*999)
813  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
814  & field_vector_dimension_type,err,error,*999)
815  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
816  & field_dp_type,err,error,*999)
817  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
818  & number_of_dimensions,err,error,*999)
819  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
820  !Constant source. Materials field components are 1 for each dimension and 1 for the constant source
821  !i.e., k and c in div(k.grad(u(x)))=c(x)
822  number_of_materials_components=number_of_dimensions+1
823  ELSE
824  !Linear source. Materials field components are 1 for each dimension and 2 for the linear source
825  !i.e., k and a and c in div(k.grad(u(x)))=a(x)u(x)+c(x)
826  number_of_materials_components=number_of_dimensions+2
827  ENDIF
828  !Set the number of materials components
829  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
830  & number_of_materials_components,err,error,*999)
831  !Default the k materials components to the geometric interpolation setup with constant interpolation
832  DO component_idx=1,number_of_dimensions
833  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
834  & component_idx,geometric_component_number,err,error,*999)
835  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
836  & component_idx,geometric_component_number,err,error,*999)
837  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
838  & component_idx,field_constant_interpolation,err,error,*999)
839  ENDDO !component_idx
840  !Default the source materials components to the first component geometric interpolation with constant interpolation
841  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
842  & 1,geometric_component_number,err,error,*999)
843  DO component_idx=number_of_dimensions+1,number_of_materials_components
844  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
845  & component_idx,geometric_component_number,err,error,*999)
846  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
847  & component_idx,field_constant_interpolation,err,error,*999)
848  ENDDO !components_idx
849  !Default the field scaling to that of the geometric field
850  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
851  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
852  ELSE
853  !Check the user specified field
854  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
855  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
856  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
857  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
858  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
859  & err,error,*999)
860  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
861  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
862  & number_of_dimensions,err,error,*999)
863  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
864  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+1, &
865  & err,error,*999)
866  ELSE
867  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+2, &
868  & err,error,*999)
869  ENDIF
870  ENDIF
871  ELSE
872  CALL flagerror("Equations set materials is not associated.",err,error,*999)
873  ENDIF
874 
876  equations_materials=>equations_set%MATERIALS
877  IF(ASSOCIATED(equations_materials)) THEN
878  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
879  !Finish creating the materials field
880  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
881  !Set the default values for the materials field
882  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
883  & number_of_dimensions,err,error,*999)
884  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
885  !Constant source. Materials field components are 1 for each dimension and 1 for the constant source
886  !i.e., k and c in div(k.grad(u(x)))=c(x)
887  number_of_materials_components=number_of_dimensions+1
888  ELSE
889  !Linear source. Materials field components are 1 for each dimension and 2 for the linear source
890  !i.e., k and a and c in div(k.grad(u(x)))=a(x)u(x)+c(x)
891  number_of_materials_components=number_of_dimensions+2
892  ENDIF
893  !First set the k values to 1.0
894  DO component_idx=1,number_of_dimensions
895  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
896  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
897  ENDDO !component_idx
898  !Now set the source values to 1.0
899  DO component_idx=number_of_dimensions+1,number_of_materials_components
900  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
901  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
902  ENDDO !component_idx
903  ENDIF
904  ELSE
905  CALL flagerror("Equations set materials is not associated.",err,error,*999)
906  ENDIF
907  CASE DEFAULT
908  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
909  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
910  & " is invalid for a velocity source Poisson equation."
911  CALL flagerror(local_error,err,error,*999)
912  END SELECT
914  SELECT CASE(equations_set%SPECIFICATION(3))
917  SELECT CASE(equations_set_setup%ACTION_TYPE)
918  !Set start action
920  IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
921  !Create the auto created source field
922  !start field creation with name 'SOURCE_FIELD'
923  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
924  & equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
925  !start creation of a new field
926  CALL field_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_general_type,err,error,*999)
927  !label the field
928  CALL field_label_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,"Source Field",err,error, &
929  & *999)
930  !define new created field to be source
931  CALL field_dependent_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
932  & field_independent_type,err,error,*999)
933  !look for decomposition rule already defined
934  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
935  & err,error,*999)
936  !apply decomposition rule found on new created field
937  CALL field_mesh_decomposition_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
938  & geometric_decomposition,err,error,*999)
939  !point new field to geometric field
940  CALL field_geometric_field_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,equations_set% &
941  & geometry%GEOMETRIC_FIELD,err,error,*999)
942  !set number of variables to 1 (1 for U)
943  source_field_number_of_variables=1
944  CALL field_number_of_variables_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
945  & source_field_number_of_variables,err,error,*999)
946  CALL field_variable_types_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
947  & (/field_u_variable_type/),err,error,*999)
948  CALL field_dimension_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
949  & field_vector_dimension_type,err,error,*999)
950  CALL field_data_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
951  & field_dp_type,err,error,*999)
952  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
953  & number_of_dimensions,err,error,*999)
954  !calculate number of components with one component for each dimension
955  source_field_number_of_components=number_of_dimensions
956  CALL field_number_of_components_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
957  & field_u_variable_type,source_field_number_of_components,err,error,*999)
958  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
959  & 1,geometric_mesh_component,err,error,*999)
960  !Default to the geometric interpolation setup
961  DO i=1,source_field_number_of_components
962  CALL field_component_mesh_component_set(equations_set%SOURCE%SOURCE_FIELD, &
963  & field_u_variable_type,i,geometric_mesh_component,err,error,*999)
964  END DO
965 
966  SELECT CASE(equations_set%SOLUTION_METHOD)
967  !Specify fem solution method
969  DO i=1,source_field_number_of_components
970  CALL field_component_interpolation_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
971  & field_u_variable_type,i,field_node_based_interpolation,err,error,*999)
972  END DO
973  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
974  & err,error,*999)
975  CALL field_scaling_type_set(equations_set%SOURCE%SOURCE_FIELD,geometric_scaling_type, &
976  & err,error,*999)
977  !Other solutions not defined yet
978  CASE DEFAULT
979  local_error="The solution method of " &
980  & //trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// " is invalid."
981  CALL flagerror(local_error,err,error,*999)
982  END SELECT
983  ELSE
984  !Check the user specified field
985  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
986  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
987  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
988  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
989  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
990  & err,error,*999)
991  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
992  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
993  & number_of_dimensions,err,error,*999)
994  !calculate number of components with one component for each dimension and one for pressure
995  source_field_number_of_components=number_of_dimensions
996  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
997  & source_field_number_of_components,err,error,*999)
998  SELECT CASE(equations_set%SOLUTION_METHOD)
1000  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
1001  & field_node_based_interpolation,err,error,*999)
1002  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
1003  & field_node_based_interpolation,err,error,*999)
1004  CASE DEFAULT
1005  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
1006  &"*",err,error))//" is invalid."
1007  CALL flagerror(local_error,err,error,*999)
1008  END SELECT
1009  ENDIF
1010  !Specify finish action
1012  IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
1013  CALL field_create_finish(equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
1014  !These 2 parameter sets will contain the fitted hermite/lagrange velocity field
1015  CALL field_parameter_set_create(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1016  & field_input_data1_set_type,err,error,*999)
1017  CALL field_parameter_set_create(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1018  & field_input_data2_set_type,err,error,*999)
1019 ! CALL FIELD_PARAMETER_SET_CREATE(EQUATIONS_SET%SOURCE%SOURCE_FIELD,FIELD_U_VARIABLE_TYPE, &
1020 ! & FIELD_INPUT_DATA3_SET_TYPE,ERR,ERROR,*999)
1021 ! CALL FIELD_PARAMETER_SET_CREATE(EQUATIONS_SET%SOURCE%SOURCE_FIELD,FIELD_U_VARIABLE_TYPE, &
1022 ! & FIELD_BOUNDARY_SET_TYPE,ERR,ERROR,*999)
1023  ENDIF
1024  CASE DEFAULT
1025  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1026  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1027  & " is invalid for a standard Navier-Poisson fluid"
1028  CALL flagerror(local_error,err,error,*999)
1029  END SELECT
1030  CASE DEFAULT
1031  local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1032  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1033  & " is invalid for a Navier-Poisson equation."
1034  CALL flagerror(local_error,err,error,*999)
1035  END SELECT
1036  !define an independent field for ALE information
1038  SELECT CASE(equations_set%SPECIFICATION(3))
1040  !do nothing!
1042  !do nothing!
1044  SELECT CASE(equations_set_setup%ACTION_TYPE)
1045  !Set start action
1047  IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
1048  !Create the auto created independent field
1049  !start field creation with name 'INDEPENDENT_FIELD'
1050  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
1051  & equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
1052  !start creation of a new field
1053  CALL field_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_general_type,err,error,*999)
1054  !label the field
1055  CALL field_label_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,"Independent Field",err,error,*999)
1056  !define new created field to be independent
1057  CALL field_dependent_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
1058  & field_independent_type,err,error,*999)
1059  !look for decomposition rule already defined
1060  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
1061  & err,error,*999)
1062  !apply decomposition rule found on new created field
1063  CALL field_mesh_decomposition_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
1064  & geometric_decomposition,err,error,*999)
1065  !point new field to geometric field
1066  CALL field_geometric_field_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,equations_set% &
1067  & geometry%GEOMETRIC_FIELD,err,error,*999)
1068  !set number of variables to 1 (1 for U)
1069  CALL field_number_of_variables_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
1070  & 1,err,error,*999)
1071  CALL field_variable_types_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
1072  & (/field_u_variable_type/),err,error,*999)
1073  CALL field_dimension_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
1074  & field_vector_dimension_type,err,error,*999)
1075  CALL field_data_type_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
1076  & field_dp_type,err,error,*999)
1077  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1078  & number_of_dimensions,err,error,*999)
1079  !calculate number of components with one component for each dimension
1080  CALL field_number_of_components_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
1081  & field_u_variable_type,number_of_dimensions,err,error,*999)
1082  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1083  & 1,geometric_mesh_component,err,error,*999)
1084  !Default to the geometric interpolation setup
1085  DO i=1,number_of_dimensions
1086  CALL field_component_mesh_component_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
1087  & field_u_variable_type,i,geometric_mesh_component,err,error,*999)
1088  END DO
1089  SELECT CASE(equations_set%SOLUTION_METHOD)
1090  !Specify fem solution method
1092  DO i=1,number_of_dimensions
1093  CALL field_component_interpolation_set_and_lock(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
1094  & field_u_variable_type,i,field_node_based_interpolation,err,error,*999)
1095  END DO
1096  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
1097  & err,error,*999)
1098  CALL field_scaling_type_set(equations_set%INDEPENDENT%INDEPENDENT_FIELD,geometric_scaling_type, &
1099  & err,error,*999)
1100  !Other solutions not defined yet
1101  CASE DEFAULT
1102  local_error="The solution method of " &
1103  & //trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// " is invalid."
1104  CALL flagerror(local_error,err,error,*999)
1105  END SELECT
1106  ELSE
1107  !Check the user specified field
1108  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1109  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1110  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
1111  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
1112  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1113  & err,error,*999)
1114  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1115  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1116  & number_of_dimensions,err,error,*999)
1117  !calculate number of components with one component for each dimension and one for pressure
1118  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
1119  & number_of_dimensions,err,error,*999)
1120  SELECT CASE(equations_set%SOLUTION_METHOD)
1122  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
1123  & field_node_based_interpolation,err,error,*999)
1124  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
1125  & field_node_based_interpolation,err,error,*999)
1126  CASE DEFAULT
1127  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
1128  &"*",err,error))//" is invalid."
1129  CALL flagerror(local_error,err,error,*999)
1130  END SELECT
1131  ENDIF
1132  !Specify finish action
1134  IF(equations_set%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN
1135  CALL field_create_finish(equations_set%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999)
1136  CALL field_parameter_set_create(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
1137  & field_mesh_displacement_set_type,err,error,*999)
1138  CALL field_parameter_set_create(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
1139  & field_mesh_velocity_set_type,err,error,*999)
1140  CALL field_parameter_set_create(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
1141  & field_boundary_set_type,err,error,*999)
1142  ENDIF
1143  CASE DEFAULT
1144  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1145  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1146  & " is invalid for a standard PPE fluid"
1147  CALL flagerror(local_error,err,error,*999)
1148  END SELECT
1149  CASE DEFAULT
1150  local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1151  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1152  & " is invalid for a PPE equation."
1153  CALL flagerror(local_error,err,error,*999)
1154  END SELECT
1156  SELECT CASE(equations_set_setup%ACTION_TYPE)
1158  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
1159  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
1160  IF(ASSOCIATED(dependent_field)) THEN
1161  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
1162  IF(ASSOCIATED(geometric_field)) THEN
1163  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
1164  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
1165 ! CASE(EQUATIONS_SET_STOKES_EQUATION_TWO_DIM_1)
1166 ! !Set analtyic function type
1167 ! EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET_STOKES_EQUATION_TWO_DIM_1
1168 ! CASE(EQUATIONS_SET_STOKES_EQUATION_TWO_DIM_2)
1169 ! !Set analtyic function type
1170 ! EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET_STOKES_EQUATION_TWO_DIM_2
1171 ! CASE(EQUATIONS_SET_STOKES_EQUATION_TWO_DIM_3)
1172 ! !Set analtyic function type
1173 ! EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET_STOKES_EQUATION_TWO_DIM_3
1174 ! CASE(EQUATIONS_SET_STOKES_EQUATION_THREE_DIM_1)
1175 ! !Set analtyic function type
1176 ! EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET_STOKES_EQUATION_THREE_DIM_1
1177 ! CASE(EQUATIONS_SET_STOKES_EQUATION_THREE_DIM_2)
1178 ! !Set analtyic function type
1179 ! EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET_STOKES_EQUATION_THREE_DIM_2
1180 ! CASE(EQUATIONS_SET_POISSON_EQUATION_THREE_DIM_3)
1181 ! !Set analtyic function type
1182 ! EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET_POISSON_EQUATION_THREE_DIM_3
1184  !Set analtyic function type
1185  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_pressure_poisson_three_dim_1
1187  !Set analtyic function type
1188  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_pressure_poisson_three_dim_2
1189  CASE DEFAULT
1190  local_error="The specified analytic function type of "// &
1191  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1192  & " is invalid for a velocity source Poisson equation."
1193  CALL flagerror(local_error,err,error,*999)
1194  END SELECT
1195  ELSE
1196  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
1197  ENDIF
1198  ELSE
1199  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
1200  ENDIF
1201  ELSE
1202  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
1203  ENDIF
1205  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
1206  analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
1207  IF(ASSOCIATED(analytic_field)) THEN
1208  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
1209  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
1210  ENDIF
1211  ENDIF
1212  ELSE
1213  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
1214  ENDIF
1215  CASE DEFAULT
1216  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1217  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1218  & " is invalid for a velocity source Poisson equation."
1219  CALL flagerror(local_error,err,error,*999)
1220  END SELECT
1222  SELECT CASE(equations_set_setup%ACTION_TYPE)
1224  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
1225  !Create the equations
1226  CALL equations_create_start(equations_set,equations,err,error,*999)
1227  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
1228  CALL equations_time_dependence_type_set(equations,equations_quasistatic,err,error,*999)
1229  ELSE
1230  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
1231  ENDIF
1233  SELECT CASE(equations_set%SOLUTION_METHOD)
1235  !Finish the creation of the equations
1236  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
1237  CALL equations_create_finish(equations,err,error,*999)
1238  !Create the equations mapping.
1239  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
1240  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
1241  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,(/field_u_variable_type/), &
1242  & err,error,*999)
1243  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
1244  CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
1245  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
1246  !Create the equations matrices
1247  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
1248  SELECT CASE(equations%SPARSITY_TYPE)
1251  & err,error,*999)
1254  & err,error,*999)
1256  & err,error,*999)
1257  CASE DEFAULT
1258  local_error="The equations matrices sparsity type of "// &
1259  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
1260  CALL flagerror(local_error,err,error,*999)
1261  END SELECT
1262  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
1264  CALL flagerror("Not implemented.",err,error,*999)
1266  CALL flagerror("Not implemented.",err,error,*999)
1268  CALL flagerror("Not implemented.",err,error,*999)
1270  CALL flagerror("Not implemented.",err,error,*999)
1272  CALL flagerror("Not implemented.",err,error,*999)
1273  CASE DEFAULT
1274  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1275  & " is invalid."
1276  CALL flagerror(local_error,err,error,*999)
1277  END SELECT
1278  CASE DEFAULT
1279  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1280  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1281  & " is invalid for a velocity source Poisson equation."
1282  CALL flagerror(local_error,err,error,*999)
1283  END SELECT
1284  CASE DEFAULT
1285  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1286  & " is invalid for a velocity source Poisson equation."
1287  CALL flagerror(local_error,err,error,*999)
1288  END SELECT
1289  ELSE
1290  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1291  & " is not a velocity source Poisson equation subtype."
1292  CALL flagerror(local_error,err,error,*999)
1293  ENDIF
1294  ELSE
1295  CALL flagerror("Equations set is not associated.",err,error,*999)
1296  ENDIF
1297 
1298  exits("Poisson_EquationsSetPressurePoissonSetup")
1299  RETURN
1300 999 errors("Poisson_EquationsSetPressurePoissonSetup",err,error)
1301  exits("Poisson_EquationsSetPressurePoissonSetup")
1302  RETURN 1
1303 
1305 
1306  !
1307  !================================================================================================================================
1308  !
1310 
1311  SUBROUTINE poisson_pre_solve_update_ppe_mesh(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
1313  !Argument variables
1314  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
1315  TYPE(solver_type), POINTER :: SOLVER
1316  INTEGER(INTG), INTENT(OUT) :: ERR
1317  TYPE(varying_string), INTENT(OUT) :: ERROR
1318  !Local Variables
1319  TYPE(solver_type), POINTER :: SOLVER_ALE_PPE, SOLVER_LAPLACE
1320  TYPE(field_type), POINTER :: INDEPENDENT_FIELD_ALE_PPE
1321  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS_ALE_PPE
1322  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING_ALE_PPE
1323  TYPE(equations_set_type), POINTER :: EQUATIONS_SET_ALE_PPE
1324  TYPE(equations_type), POINTER :: EQUATIONS
1325  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
1326  TYPE(varying_string) :: LOCAL_ERROR
1327  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
1328  TYPE(domain_type), POINTER :: DOMAIN
1329  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
1330 
1331  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT,ALPHA
1332  REAL(DP), POINTER :: MESH_DISPLACEMENT_VALUES(:)
1333  INTEGER(INTG) :: NUMBER_OF_DIMENSIONS_ALE_PPE,GEOMETRIC_MESH_COMPONENT
1334  INTEGER(INTG) :: INPUT_TYPE,INPUT_OPTION,component_idx,deriv_idx,local_ny,node_idx,variable_idx,variable_type
1335 
1336  enters("POISSON_PRE_SOLVE_UPDATE_PPE_MESH",err,error,*999)
1337 
1338  IF(ASSOCIATED(control_loop)) THEN
1339  CALL control_loop_current_times_get(control_loop%PARENT_LOOP,current_time,time_increment,err,error,*999)
1340  NULLIFY(solver_laplace)
1341  NULLIFY(solver_ale_ppe)
1342  IF(ASSOCIATED(solver)) THEN
1343  IF(ASSOCIATED(control_loop%PARENT_LOOP%PROBLEM)) THEN
1344  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
1345  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1346  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
1347  CALL flagerror("Problem specification must have three entries for a Poisson problem.",err,error,*999)
1348  END IF
1349  SELECT CASE(control_loop%PARENT_LOOP%PROBLEM%SPECIFICATION(3))
1351  ! do nothing ???
1353  ! do nothing ???
1355  !Update mesh within the dynamic solver
1356  IF(solver%SOLVE_TYPE==solver_linear_type) THEN
1357  !Get the independent field for the ALE PPE problem
1358  CALL solvers_solver_get(solver%SOLVERS,1,solver_ale_ppe,err,error,*999)
1359  solver_equations_ale_ppe=>solver_ale_ppe%SOLVER_EQUATIONS
1360  IF(ASSOCIATED(solver_equations_ale_ppe)) THEN
1361  solver_mapping_ale_ppe=>solver_equations_ale_ppe%SOLVER_MAPPING
1362  IF(ASSOCIATED(solver_mapping_ale_ppe)) THEN
1363  equations_set_ale_ppe=>solver_mapping_ale_ppe%EQUATIONS_SETS(1)%PTR
1364  IF(ASSOCIATED(equations_set_ale_ppe)) THEN
1365  independent_field_ale_ppe=>equations_set_ale_ppe%INDEPENDENT%INDEPENDENT_FIELD
1366  ELSE
1367  CALL flagerror("ALE PPE equations set is not associated.",err,error,*999)
1368  END IF
1369  !Get the data
1370  CALL field_number_of_components_get(equations_set_ale_ppe%GEOMETRY%GEOMETRIC_FIELD, &
1371  & field_u_variable_type,number_of_dimensions_ale_ppe,err,error,*999)
1372 !\todo: Introduce flags set by the user (42/1 only for testings purpose)
1373  !Copy input to PPE' independent field
1374  input_type=42
1375  input_option=1
1376  NULLIFY(mesh_displacement_values)
1377  CALL field_parameter_set_data_get(equations_set_ale_ppe%INDEPENDENT%INDEPENDENT_FIELD, &
1378  & field_u_variable_type,field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
1379  CALL fluid_mechanics_io_read_data(solver_linear_type,mesh_displacement_values, &
1380  & number_of_dimensions_ale_ppe,input_type,input_option,control_loop%PARENT_LOOP%TIME_LOOP% &
1381  & iteration_number,1.0_dp)
1382  CALL field_parameter_set_update_start(equations_set_ale_ppe%INDEPENDENT%INDEPENDENT_FIELD, &
1383  & field_u_variable_type,field_mesh_displacement_set_type,err,error,*999)
1384  CALL field_parameter_set_update_finish(equations_set_ale_ppe%INDEPENDENT%INDEPENDENT_FIELD, &
1385  & field_u_variable_type,field_mesh_displacement_set_type,err,error,*999)
1386  ELSE
1387  CALL flagerror("ALE PPE solver mapping is not associated.",err,error,*999)
1388  END IF
1389  ELSE
1390  CALL flagerror("ALE PPE solver equations are not associated.",err,error,*999)
1391  END IF
1392  !Use calculated values to update mesh
1393  CALL field_component_mesh_component_get(equations_set_ale_ppe%GEOMETRY%GEOMETRIC_FIELD, &
1394  & field_u_variable_type,1,geometric_mesh_component,err,error,*999)
1395 ! CALL FIELD_PARAMETER_SET_DATA_GET(INDEPENDENT_FIELD_ALE_PPE,FIELD_U_VARIABLE_TYPE, &
1396 ! & FIELD_MESH_DISPLACEMENT_SET_TYPE,MESH_DISPLACEMENT_VALUES,ERR,ERROR,*999)
1397  equations=>solver_mapping_ale_ppe%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
1398  IF(ASSOCIATED(equations)) THEN
1399  equations_mapping=>equations%EQUATIONS_MAPPING
1400  IF(ASSOCIATED(equations_mapping)) THEN
1401  DO variable_idx=1,equations_set_ale_ppe%DEPENDENT%DEPENDENT_FIELD%NUMBER_OF_VARIABLES
1402  variable_type=equations_set_ale_ppe%DEPENDENT%DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
1403  field_variable=>equations_set_ale_ppe%GEOMETRY%GEOMETRIC_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
1404  IF(ASSOCIATED(field_variable)) THEN
1405  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
1406  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
1407  IF(ASSOCIATED(domain)) THEN
1408  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
1409  domain_nodes=>domain%TOPOLOGY%NODES
1410  IF(ASSOCIATED(domain_nodes)) THEN
1411  !Loop over the local nodes excluding the ghosts.
1412  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
1413  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
1414  !Default to version 1 of each node derivative
1415  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
1416  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
1417  CALL field_parameter_set_add_local_dof(equations_set_ale_ppe%GEOMETRY%GEOMETRIC_FIELD, &
1418  & field_u_variable_type,field_values_set_type,local_ny, &
1419  & mesh_displacement_values(local_ny),err,error,*999)
1420  ENDDO !deriv_idx
1421  ENDDO !node_idx
1422  ENDIF
1423  ENDIF
1424  ENDIF
1425  ENDDO !component_idx
1426  ENDIF
1427  ENDDO !variable_idx
1428  ELSE
1429  CALL flagerror("Equations mapping is not associated.",err,error,*999)
1430  END IF
1431  ELSE
1432  CALL flagerror("Equations are not associated.",err,error,*999)
1433  END IF
1434  CALL field_parameter_set_update_start(equations_set_ale_ppe%GEOMETRY%GEOMETRIC_FIELD, &
1435  & field_u_variable_type,field_values_set_type,err,error,*999)
1436  CALL field_parameter_set_update_finish(equations_set_ale_ppe%GEOMETRY%GEOMETRIC_FIELD, &
1437  & field_u_variable_type,field_values_set_type,err,error,*999)
1438  !Now use displacement values to calculate velocity values
1439  time_increment=control_loop%PARENT_LOOP%TIME_LOOP%TIME_INCREMENT
1440  alpha=1.0_dp/time_increment
1441  CALL field_parameter_sets_copy(independent_field_ale_ppe,field_u_variable_type, &
1442  & field_mesh_displacement_set_type,field_mesh_velocity_set_type,alpha,err,error,*999)
1443  ELSE
1444  CALL flagerror("Mesh motion calculation not successful for ALE problem.",err,error,*999)
1445  END IF
1446  CASE DEFAULT
1447  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
1448  & " is not valid for a PPE equation fluid type of a fluid mechanics problem class."
1449  CALL flagerror(local_error,err,error,*999)
1450  END SELECT
1451  ELSE
1452  CALL flagerror("Problem is not associated.",err,error,*999)
1453  ENDIF
1454  ELSE
1455  CALL flagerror("Solver is not associated.",err,error,*999)
1456  ENDIF
1457  ELSE
1458  CALL flagerror("Control loop is not associated.",err,error,*999)
1459  ENDIF
1460  exits("POISSON_PRE_SOLVE_UPDATE_PPE_MESH")
1461  RETURN
1462 999 errorsexits("POISSON_PRE_SOLVE_UPDATE_PPE_MESH",err,error)
1463  RETURN 1
1464  END SUBROUTINE poisson_pre_solve_update_ppe_mesh
1465 
1466  !
1467  !================================================================================================================================
1468  !
1470 
1471  SUBROUTINE poisson_pre_solve_update_ppe_source(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
1473  !Argument variables
1474  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
1475  TYPE(solver_type), POINTER :: SOLVER
1476  INTEGER(INTG), INTENT(OUT) :: ERR
1477  TYPE(varying_string), INTENT(OUT) :: ERROR
1478  !Local Variables
1479  TYPE(solver_type), POINTER :: SOLVER_FITTED, SOLVER_PPE
1480  TYPE(field_type), POINTER :: DEPENDENT_FIELD_FITTED,SOURCE_FIELD_PPE
1481  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS_FITTED, SOLVER_EQUATIONS_PPE
1482  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING_FITTED, SOLVER_MAPPING_PPE
1483  TYPE(equations_set_type), POINTER :: EQUATIONS_SET_FITTED, EQUATIONS_SET_PPE
1484 ! TYPE(EQUATIONS_TYPE), POINTER :: EQUATIONS
1485 ! TYPE(EQUATIONS_MAPPING_TYPE), POINTER :: EQUATIONS_MAPPING
1486  TYPE(varying_string) :: LOCAL_ERROR
1487 ! TYPE(FIELD_VARIABLE_TYPE), POINTER :: FIELD_VARIABLE
1488 ! TYPE(DOMAIN_TYPE), POINTER :: DOMAIN
1489 ! TYPE(DOMAIN_NODES_TYPE), POINTER :: DOMAIN_NODES
1490 
1491  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
1492  INTEGER(INTG) :: I,NUMBER_OF_DIMENSIONS_PPE,NUMBER_OF_DIMENSIONS_FITTED,GEOMETRIC_MESH_COMPONENT
1493 
1494  enters("POISSON_PRE_SOLVE_UPDATE_PPE_SOURCE",err,error,*999)
1495 
1496  IF(ASSOCIATED(control_loop)) THEN
1497  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
1498  NULLIFY(solver_ppe)
1499  NULLIFY(solver_fitted)
1500  IF(ASSOCIATED(solver)) THEN
1501  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
1502  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
1503  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1504  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
1505  CALL flagerror("Problem specification must have three entries for a Poisson problem.",err,error,*999)
1506  END IF
1507  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
1509  ! do nothing ???
1511  ! do nothing ???
1513  ! do nothing ???
1515  IF(solver%GLOBAL_NUMBER==1) THEN
1516  !do nothing for the fitting process
1517  ELSE IF(solver%GLOBAL_NUMBER==2) THEN
1518  !Get the dependent field for the three component Laplace problem
1519  CALL solvers_solver_get(solver%SOLVERS,1,solver_fitted,err,error,*999)
1520  solver_equations_fitted=>solver_fitted%SOLVER_EQUATIONS
1521  IF(ASSOCIATED(solver_equations_fitted)) THEN
1522  solver_mapping_fitted=>solver_equations_fitted%SOLVER_MAPPING
1523  IF(ASSOCIATED(solver_mapping_fitted)) THEN
1524  equations_set_fitted=>solver_mapping_fitted%EQUATIONS_SETS(1)%PTR
1525  IF(ASSOCIATED(equations_set_fitted)) THEN
1526  dependent_field_fitted=>equations_set_fitted%DEPENDENT%DEPENDENT_FIELD
1527  ELSE
1528  CALL flagerror("Fitted equations set is not associated.",err,error,*999)
1529  END IF
1530  CALL field_number_of_components_get(equations_set_fitted%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1531  & number_of_dimensions_fitted,err,error,*999)
1532  ELSE
1533  CALL flagerror("Fitted solver mapping is not associated.",err,error,*999)
1534  END IF
1535  ELSE
1536  CALL flagerror("Fitted solver equations are not associated.",err,error,*999)
1537  END IF
1538  !Get the source field for the PPE problem
1539  CALL solvers_solver_get(solver%SOLVERS,2,solver_ppe,err,error,*999)
1540  solver_equations_ppe=>solver_ppe%SOLVER_EQUATIONS
1541  IF(ASSOCIATED(solver_equations_ppe)) THEN
1542  solver_mapping_ppe=>solver_equations_ppe%SOLVER_MAPPING
1543  IF(ASSOCIATED(solver_mapping_ppe)) THEN
1544  equations_set_ppe=>solver_mapping_ppe%EQUATIONS_SETS(1)%PTR
1545  IF(ASSOCIATED(equations_set_ppe)) THEN
1546  source_field_ppe=>equations_set_ppe%SOURCE%SOURCE_FIELD
1547  ELSE
1548  CALL flagerror("PPE equations set is not associated.",err,error,*999)
1549  END IF
1550  CALL field_number_of_components_get(equations_set_ppe%GEOMETRY%GEOMETRIC_FIELD, &
1551  & field_u_variable_type,number_of_dimensions_ppe,err,error,*999)
1552  ELSE
1553  CALL flagerror("PPE solver mapping is not associated.",err,error,*999)
1554  END IF
1555  ELSE
1556  CALL flagerror("PPE solver equations are not associated.",err,error,*999)
1557  END IF
1558  !Copy result from FITTING to PPE's source field
1559  IF(number_of_dimensions_ppe==number_of_dimensions_fitted) THEN
1560  IF(control_loop%TIME_LOOP%ITERATION_NUMBER/=1) THEN
1561  DO i=1,number_of_dimensions_ppe
1562  CALL write_string(general_output_type,"Save old source field... ",err,error,*999)
1563  CALL field_parameterstofieldparameterscopy(source_field_ppe, &
1564  & field_u_variable_type,field_input_data1_set_type,i,source_field_ppe, &
1565  & field_u_variable_type,field_input_data2_set_type,i,err,error,*999)
1566  END DO
1567  ENDIF
1568  CALL write_string(general_output_type,"Update new source field... ",err,error,*999)
1569  DO i=1,number_of_dimensions_ppe
1570  CALL field_parameterstofieldparameterscopy(dependent_field_fitted, &
1571  & field_u_variable_type,field_values_set_type,i,source_field_ppe, &
1572  & field_u_variable_type,field_input_data1_set_type,i,err,error,*999)
1573  END DO
1574  IF(control_loop%TIME_LOOP%ITERATION_NUMBER==1) THEN
1575  CALL write_string(general_output_type,"Old source field is new source field!... ",err,error,*999)
1576  DO i=1,number_of_dimensions_ppe
1577  CALL field_parameterstofieldparameterscopy(source_field_ppe, &
1578  & field_u_variable_type,field_input_data1_set_type,i,source_field_ppe, &
1579  & field_u_variable_type,field_input_data2_set_type,i,err,error,*999)
1580  END DO
1581  ENDIF
1582  ELSE
1583  CALL flagerror("Dimension of FITTED and PPE equations set is not consistent.",err,error,*999)
1584  END IF
1585  !Use calculated values to update mesh
1586  CALL field_component_mesh_component_get(equations_set_ppe%GEOMETRY%GEOMETRIC_FIELD, &
1587  & field_u_variable_type,1,geometric_mesh_component,err,error,*999)
1588  CALL field_parameter_set_update_start(equations_set_ppe%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1589  & field_input_data1_set_type,err,error,*999)
1590  CALL field_parameter_set_update_finish(equations_set_ppe%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1591  & field_input_data1_set_type,err,error,*999)
1592  CALL field_parameter_set_update_start(equations_set_ppe%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1593  & field_input_data2_set_type,err,error,*999)
1594  CALL field_parameter_set_update_finish(equations_set_ppe%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1595  & field_input_data2_set_type,err,error,*999)
1596  ELSE
1597  CALL flagerror("Mesh update is not defined for non-dynamic problems.",err,error,*999)
1598  END IF
1599  CASE DEFAULT
1600  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
1601  & " is not valid for a PPE equation fluid type of a fluid mechanics problem class."
1602  CALL flagerror(local_error,err,error,*999)
1603  END SELECT
1604  ELSE
1605  CALL flagerror("Problem is not associated.",err,error,*999)
1606  ENDIF
1607  ELSE
1608  CALL flagerror("Solver is not associated.",err,error,*999)
1609  ENDIF
1610  ELSE
1611  CALL flagerror("Control loop is not associated.",err,error,*999)
1612  ENDIF
1613  exits("POISSON_PRE_SOLVE_UPDATE_PPE_SOURCE")
1614  RETURN
1615 999 errorsexits("POISSON_PRE_SOLVE_UPDATE_PPE_SOURCE",err,error)
1616  RETURN 1
1618 
1619  !
1620  !================================================================================================================================
1621  !
1622 
1623 
1625  SUBROUTINE poisson_equationssetlinearsourcesetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
1627  !Argument variables
1628  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1629  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
1630  INTEGER(INTG), INTENT(OUT) :: ERR
1631  TYPE(varying_string), INTENT(OUT) :: ERROR
1632  !Local Variables
1633  INTEGER(INTG) :: component_idx,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE,NUMBER_OF_DIMENSIONS, &
1634  & NUMBER_OF_MATERIALS_COMPONENTS, GEOMETRIC_MESH_COMPONENT,SOURCE_FIELD_NUMBER_OF_COMPONENTS,I, &
1635  & SOURCE_FIELD_NUMBER_OF_VARIABLES
1636  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
1637  TYPE(equations_type), POINTER :: EQUATIONS
1638  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
1639  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
1640  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
1641  TYPE(field_type), POINTER :: analytic_field,DEPENDENT_FIELD,geometric_field
1642  TYPE(varying_string) :: LOCAL_ERROR
1643 
1644  enters("POISSON_EQUATION_EQUATION_SET_LINEAR_SOURCE_SETUP",err,error,*999)
1645 
1646  NULLIFY(equations)
1647  NULLIFY(equations_mapping)
1648  NULLIFY(equations_matrices)
1649  NULLIFY(geometric_decomposition)
1650 
1651  IF(ASSOCIATED(equations_set)) THEN
1652  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
1653  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
1654  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
1655  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
1656  & err,error,*999)
1657  END IF
1658  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype.OR. &
1659  & equations_set%SPECIFICATION(3)==equations_set_linear_source_poisson_subtype) THEN
1660  SELECT CASE(equations_set_setup%SETUP_TYPE)
1662  SELECT CASE(equations_set_setup%ACTION_TYPE)
1666 !!TODO: Check valid setup
1667  CASE DEFAULT
1668  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1669  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1670  & " is invalid for a linear source Poisson equation."
1671  CALL flagerror(local_error,err,error,*999)
1672  END SELECT
1674  !Do nothing???
1676  SELECT CASE(equations_set_setup%ACTION_TYPE)
1678  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
1679  !Create the auto created dependent field
1680  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
1681  & dependent_field,err,error,*999)
1682  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
1683  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
1684  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1685  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
1686  & err,error,*999)
1687  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
1688  & geometric_field,err,error,*999)
1689  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
1690  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,(/field_u_variable_type, &
1691  & field_deludeln_variable_type/),err,error,*999)
1692  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
1693  & field_scalar_dimension_type,err,error,*999)
1694  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
1695  & field_scalar_dimension_type,err,error,*999)
1696  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
1697  & field_dp_type,err,error,*999)
1698  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
1699  & field_dp_type,err,error,*999)
1700  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
1701  & err,error,*999)
1702  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
1703  & err,error,*999)
1704  !Default to the geometric interpolation setup
1705  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
1706  & geometric_component_number,err,error,*999)
1707  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
1708  & geometric_component_number,err,error,*999)
1709  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
1710  & geometric_component_number,err,error,*999)
1711  SELECT CASE(equations_set%SOLUTION_METHOD)
1713  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
1714  & field_node_based_interpolation,err,error,*999)
1715  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
1716  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
1717  !Default the scaling to the geometric field scaling
1718  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1719  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
1721  CALL flagerror("Not implemented.",err,error,*999)
1723  CALL flagerror("Not implemented.",err,error,*999)
1725  CALL flagerror("Not implemented.",err,error,*999)
1727  CALL flagerror("Not implemented.",err,error,*999)
1729  CALL flagerror("Not implemented.",err,error,*999)
1730  CASE DEFAULT
1731  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1732  & " is invalid."
1733  CALL flagerror(local_error,err,error,*999)
1734  END SELECT
1735  ELSE
1736  !Check the user specified field
1737  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1738  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
1739  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
1740  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type,field_deludeln_variable_type/), &
1741  & err,error,*999)
1742  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type,err,error,*999)
1743  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
1744  & err,error,*999)
1745  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1746  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
1747  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
1748  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1,err,error,*999)
1749  SELECT CASE(equations_set%SOLUTION_METHOD)
1751  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
1752  & field_node_based_interpolation,err,error,*999)
1753  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
1754  & field_node_based_interpolation,err,error,*999)
1756  CALL flagerror("Not implemented.",err,error,*999)
1758  CALL flagerror("Not implemented.",err,error,*999)
1760  CALL flagerror("Not implemented.",err,error,*999)
1762  CALL flagerror("Not implemented.",err,error,*999)
1764  CALL flagerror("Not implemented.",err,error,*999)
1765  CASE DEFAULT
1766  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1767  & " is invalid."
1768  CALL flagerror(local_error,err,error,*999)
1769  END SELECT
1770  ENDIF
1772  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
1773  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
1774  ENDIF
1775  CASE DEFAULT
1776  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1777  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1778  & " is invalid for a linear source Poisson equation"
1779  CALL flagerror(local_error,err,error,*999)
1780  END SELECT
1782  SELECT CASE(equations_set_setup%ACTION_TYPE)
1784  equations_materials=>equations_set%MATERIALS
1785  IF(ASSOCIATED(equations_materials)) THEN
1786  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
1787  !Create the auto created materials field
1788  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
1789  & materials_field,err,error,*999)
1790  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
1791  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
1792  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1793  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
1794  & err,error,*999)
1795  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
1796  & geometric_field,err,error,*999)
1797  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
1798  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,(/field_u_variable_type/), &
1799  & err,error,*999)
1800  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1801  & field_vector_dimension_type,err,error,*999)
1802  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1803  & field_dp_type,err,error,*999)
1804  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1805  & number_of_dimensions,err,error,*999)
1806  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
1807  !Constant source. Materials field components are 1 for each dimension and 1 for the constant source
1808  !i.e., k and c in div(k.grad(u(x)))=c(x)
1809  number_of_materials_components=number_of_dimensions+1
1810  ELSE
1811  !Linear source. Materials field components are 1 for each dimension and 2 for the linear source
1812  !i.e., k and a and c in div(k.grad(u(x)))=a(x)u(x)+c(x)
1813  number_of_materials_components=number_of_dimensions+2
1814  ENDIF
1815  !Set the number of materials components
1816  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1817  & number_of_materials_components,err,error,*999)
1818  !Default the k materials components to the geometric interpolation setup with constant interpolation
1819  DO component_idx=1,number_of_dimensions
1820  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1821  & component_idx,geometric_component_number,err,error,*999)
1822  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1823  & component_idx,geometric_component_number,err,error,*999)
1824  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1825  & component_idx,field_constant_interpolation,err,error,*999)
1826  ENDDO !component_idx
1827  !Default the source materials components to the first component geometric interpolation with constant interpolation
1828  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1829  & 1,geometric_component_number,err,error,*999)
1830  DO component_idx=number_of_dimensions+1,number_of_materials_components
1831  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1832  & component_idx,geometric_component_number,err,error,*999)
1833  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1834  & component_idx,field_constant_interpolation,err,error,*999)
1835  ENDDO !components_idx
1836  !Default the field scaling to that of the geometric field
1837  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1838  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
1839  ELSE
1840  !Check the user specified field
1841  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
1842  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1843  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
1844  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
1845  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1846  & err,error,*999)
1847  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1848  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1849  & number_of_dimensions,err,error,*999)
1850  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
1851  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+1, &
1852  & err,error,*999)
1853  ELSE
1854  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+2, &
1855  & err,error,*999)
1856  ENDIF
1857  ENDIF
1858  ELSE
1859  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1860  ENDIF
1861 
1863  equations_materials=>equations_set%MATERIALS
1864  IF(ASSOCIATED(equations_materials)) THEN
1865  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
1866  !Finish creating the materials field
1867  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
1868  !Set the default values for the materials field
1869  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1870  & number_of_dimensions,err,error,*999)
1871  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
1872  !Constant source. Materials field components are 1 for each dimension and 1 for the constant source
1873  !i.e., k and c in div(k.grad(u(x)))=c(x)
1874  number_of_materials_components=number_of_dimensions+1
1875  ELSE
1876  !Linear source. Materials field components are 1 for each dimension and 2 for the linear source
1877  !i.e., k and a and c in div(k.grad(u(x)))=a(x)u(x)+c(x)
1878  number_of_materials_components=number_of_dimensions+2
1879  ENDIF
1880  !First set the k values to 1.0
1881  DO component_idx=1,number_of_dimensions
1882  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1883  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
1884  ENDDO !component_idx
1885  !Now set the source values to 1.0
1886  DO component_idx=number_of_dimensions+1,number_of_materials_components
1887  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1888  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
1889  ENDDO !component_idx
1890  ENDIF
1891  ELSE
1892  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1893  ENDIF
1894  CASE DEFAULT
1895  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1896  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1897  & " is invalid for a linear source Poisson equation."
1898  CALL flagerror(local_error,err,error,*999)
1899  END SELECT
1901  SELECT CASE(equations_set%SPECIFICATION(3))
1903  SELECT CASE(equations_set_setup%ACTION_TYPE)
1904  !Set start action
1906  IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
1907  !Create the auto created source field
1908  !start field creation with name 'SOURCE_FIELD'
1909  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
1910  & equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
1911  !start creation of a new field
1912  CALL field_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_general_type,err,error,*999)
1913  !label the field
1914  CALL field_label_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,"Source Field",err,error, &
1915  & *999)
1916  !define new created field to be source
1917  CALL field_dependent_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
1918  & field_independent_type,err,error,*999)
1919  !look for decomposition rule already defined
1920  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
1921  & err,error,*999)
1922  !apply decomposition rule found on new created field
1923  CALL field_mesh_decomposition_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
1924  & geometric_decomposition,err,error,*999)
1925  !point new field to geometric field
1926  CALL field_geometric_field_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,equations_set% &
1927  & geometry%GEOMETRIC_FIELD,err,error,*999)
1928  !set number of variables to 1 (1 for U)
1929  source_field_number_of_variables=1
1930  CALL field_number_of_variables_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
1931  & source_field_number_of_variables,err,error,*999)
1932  CALL field_variable_types_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
1933  & (/field_u_variable_type/),err,error,*999)
1934  CALL field_dimension_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1935  & field_vector_dimension_type,err,error,*999)
1936  CALL field_data_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
1937  & field_dp_type,err,error,*999)
1938  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1939  & number_of_dimensions,err,error,*999)
1940  !calculate number of components with one component for each dimension
1941  source_field_number_of_components=number_of_dimensions
1942  CALL field_number_of_components_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
1943  & field_u_variable_type,source_field_number_of_components,err,error,*999)
1944  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1945  & 1,geometric_mesh_component,err,error,*999)
1946  !Default to the geometric interpolation setup
1947  DO i=1,source_field_number_of_components
1948  CALL field_component_mesh_component_set(equations_set%SOURCE%SOURCE_FIELD, &
1949  & field_u_variable_type,i,geometric_mesh_component,err,error,*999)
1950  END DO
1951  SELECT CASE(equations_set%SOLUTION_METHOD)
1952  !Specify fem solution method
1954  DO i=1,source_field_number_of_components
1955  CALL field_component_interpolation_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
1956  & field_u_variable_type,i,field_node_based_interpolation,err,error,*999)
1957  END DO
1958  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
1959  & err,error,*999)
1960  CALL field_scaling_type_set(equations_set%SOURCE%SOURCE_FIELD,geometric_scaling_type, &
1961  & err,error,*999)
1962  !Other solutions not defined yet
1963  CASE DEFAULT
1964  local_error="The solution method of " &
1965  & //trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// " is invalid."
1966  CALL flagerror(local_error,err,error,*999)
1967  END SELECT
1968  ELSE
1969  !Check the user specified field
1970  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1971  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1972  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
1973  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
1974  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1975  & err,error,*999)
1976  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1977  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1978  & number_of_dimensions,err,error,*999)
1979  !calculate number of components with one component for each dimension and one for pressure
1980  source_field_number_of_components=number_of_dimensions
1981  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
1982  & source_field_number_of_components,err,error,*999)
1983  SELECT CASE(equations_set%SOLUTION_METHOD)
1985  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
1986  & field_node_based_interpolation,err,error,*999)
1987  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
1988  & field_node_based_interpolation,err,error,*999)
1989  CASE DEFAULT
1990  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
1991  &"*",err,error))//" is invalid."
1992  CALL flagerror(local_error,err,error,*999)
1993  END SELECT
1994  ENDIF
1995  !Specify finish action
1997  IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
1998  CALL field_create_finish(equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
1999  ENDIF
2000  CASE DEFAULT
2001  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2002  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2003  & " is invalid for a linear Poisson subtype"
2004  CALL flagerror(local_error,err,error,*999)
2005  END SELECT
2006  CASE DEFAULT
2007  local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2008  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2009  & " is invalid for a linear Poisson equation."
2010  CALL flagerror(local_error,err,error,*999)
2011  END SELECT
2013  SELECT CASE(equations_set_setup%ACTION_TYPE)
2015  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
2016  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
2017  IF(ASSOCIATED(dependent_field)) THEN
2018  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
2019  IF(ASSOCIATED(geometric_field)) THEN
2020  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
2021  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
2023  !Check that we have an constant source
2024  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
2025  !Check that we are in 2D
2026  IF(number_of_dimensions/=2) THEN
2027  local_error="The number of geometric dimensions of "// &
2028  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2029  & " is invalid. The analytic function type of "// &
2030  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2031  & " requires that there be 2 geometric dimensions."
2032  CALL flagerror(local_error,err,error,*999)
2033  ENDIF
2034  !Create analytic field if required
2035  !Set analtyic function type
2036  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_poisson_equation_two_dim_1
2037  ELSE
2038  local_error="The equations set subtype of "// &
2039  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2040  & " is invalid. The analytic function type of "// &
2041  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2042  & " requires that the equations set subtype be an exponential source Poisson equation."
2043  CALL flagerror(local_error,err,error,*999)
2044  ENDIF
2046  !Check that we have an exponential source
2047  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
2048  !Check that we are in 3D
2049  IF(number_of_dimensions/=3) THEN
2050  local_error="The number of geometric dimensions of "// &
2051  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2052  & " is invalid. The analytic function type of "// &
2053  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2054  & " requires that there be 3 geometric dimensions."
2055  CALL flagerror(local_error,err,error,*999)
2056  ENDIF
2057  !Create analytic field if required
2058  !Set analytic function type
2059  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_poisson_equation_three_dim_1
2060  ELSE
2061  local_error="The equations set subtype of "// &
2062  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2063  & " is invalid. The analytic function type of "// &
2064  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2065  & " requires that the equations set subtype be an exponential source Poisson equation."
2066  CALL flagerror(local_error,err,error,*999)
2067  ENDIF
2068  CASE DEFAULT
2069  local_error="The specified analytic function type of "// &
2070  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2071  & " is invalid for a linear source Poisson equation."
2072  CALL flagerror(local_error,err,error,*999)
2073  END SELECT
2074  ELSE
2075  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
2076  ENDIF
2077  ELSE
2078  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
2079  ENDIF
2080  ELSE
2081  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2082  ENDIF
2084  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
2085  analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
2086  IF(ASSOCIATED(analytic_field)) THEN
2087  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
2088  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
2089  ENDIF
2090  ENDIF
2091  ELSE
2092  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
2093  ENDIF
2094  CASE DEFAULT
2095  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2096  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2097  & " is invalid for a linear source Poisson equation."
2098  CALL flagerror(local_error,err,error,*999)
2099  END SELECT
2101  SELECT CASE(equations_set_setup%ACTION_TYPE)
2103  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
2104  !Create the equations
2105  CALL equations_create_start(equations_set,equations,err,error,*999)
2106  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
2107  CALL equations_time_dependence_type_set(equations,equations_static,err,error,*999)
2108  ELSE
2109  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2110  ENDIF
2112  SELECT CASE(equations_set%SOLUTION_METHOD)
2114  !Finish the creation of the equations
2115  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
2116  CALL equations_create_finish(equations,err,error,*999)
2117  !Create the equations mapping.
2118  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
2119  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
2120  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,(/field_u_variable_type/), &
2121  & err,error,*999)
2122  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
2123  CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
2124  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
2125  !Create the equations matrices
2126  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
2127  SELECT CASE(equations%SPARSITY_TYPE)
2130  & err,error,*999)
2133  & err,error,*999)
2135  & err,error,*999)
2136  CASE DEFAULT
2137  local_error="The equations matrices sparsity type of "// &
2138  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
2139  CALL flagerror(local_error,err,error,*999)
2140  END SELECT
2141  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
2143  CALL flagerror("Not implemented.",err,error,*999)
2145  CALL flagerror("Not implemented.",err,error,*999)
2147  CALL flagerror("Not implemented.",err,error,*999)
2149  CALL flagerror("Not implemented.",err,error,*999)
2151  CALL flagerror("Not implemented.",err,error,*999)
2152  CASE DEFAULT
2153  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2154  & " is invalid."
2155  CALL flagerror(local_error,err,error,*999)
2156  END SELECT
2157  CASE DEFAULT
2158  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2159  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2160  & " is invalid for a linear source Poisson equation."
2161  CALL flagerror(local_error,err,error,*999)
2162  END SELECT
2163  CASE DEFAULT
2164  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2165  & " is invalid for a linear source Poisson equation."
2166  CALL flagerror(local_error,err,error,*999)
2167  END SELECT
2168  ELSE
2169  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2170  & " is not a linear source Poisson equation subtype."
2171  CALL flagerror(local_error,err,error,*999)
2172  ENDIF
2173  ELSE
2174  CALL flagerror("Equations set is not associated.",err,error,*999)
2175  ENDIF
2176 
2177  exits("Poisson_EquationsSetLinearSourceSetup")
2178  RETURN
2179 999 errors("Poisson_EquationsSetLinearSourceSetup",err,error)
2180  exits("Poisson_EquationsSetLinearSourceSetup")
2181  RETURN 1
2182 
2184 
2185  !
2186  !================================================================================================================================
2187  !
2188 
2190  SUBROUTINE poisson_equationssetextracellularbidomainsetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
2192  !Argument variables
2193  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2194  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
2195  INTEGER(INTG), INTENT(OUT) :: ERR
2196  TYPE(varying_string), INTENT(OUT) :: ERROR
2197  !Local Variables
2198  INTEGER(INTG) :: component_idx,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE,NUMBER_OF_DIMENSIONS, &
2199  & NUMBER_OF_MATERIALS_COMPONENTS, GEOMETRIC_MESH_COMPONENT,SOURCE_FIELD_NUMBER_OF_COMPONENTS, &
2200  & SOURCE_FIELD_NUMBER_OF_VARIABLES
2201  TYPE(field_type), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,GEOMETRIC_FIELD
2202  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
2203  TYPE(equations_type), POINTER :: EQUATIONS
2204  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
2205  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
2206  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
2207  TYPE(varying_string) :: LOCAL_ERROR
2208 
2209  enters("Poisson_EquationsSetExtracellularBidomainSetup",err,error,*999)
2210 
2211  NULLIFY(equations)
2212  NULLIFY(equations_mapping)
2213  NULLIFY(equations_matrices)
2214  NULLIFY(geometric_decomposition)
2215 
2216  IF(ASSOCIATED(equations_set)) THEN
2217  IF(equations_set%SPECIFICATION(3)==equations_set_extracellular_bidomain_poisson_subtype) THEN
2218  SELECT CASE(equations_set_setup%SETUP_TYPE)
2220  SELECT CASE(equations_set_setup%ACTION_TYPE)
2224 !!TODO: Check valid setup
2225  CASE DEFAULT
2226  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2227  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2228  & " is invalid for an extracellular bidomain Poisson equation."
2229  CALL flagerror(local_error,err,error,*999)
2230  END SELECT
2231  !GEOMETRY
2233  !Do nothing???
2234  !DEPENDENT
2236  SELECT CASE(equations_set_setup%ACTION_TYPE)
2238  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
2239  !Create the auto created dependent field
2240  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
2241  & dependent_field,err,error,*999)
2242  CALL field_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,"Dependent Field",err,error,*999)
2243  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
2244  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
2245  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
2246  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
2247  & err,error,*999)
2248  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
2249  & geometric_field,err,error,*999)
2250  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
2251  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,(/field_u_variable_type, &
2252  & field_deludeln_variable_type/),err,error,*999)
2253 
2254  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,"Phi",err,error,*999)
2255  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,"del Phi/del n", &
2256  & err,error,*999)
2257 
2258  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2259  & field_scalar_dimension_type,err,error,*999)
2260  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
2261  & field_scalar_dimension_type,err,error,*999)
2262 
2263  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2264  & field_dp_type,err,error,*999)
2265  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
2266  & field_dp_type,err,error,*999)
2267 
2268  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2269  & err,error,*999)
2270  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
2271  & err,error,*999)
2272 
2273  CALL field_component_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1,"Phi",err,error,*999)
2274  CALL field_component_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
2275  & "del Phi/del n",err,error,*999)
2276 
2277  !Default to the geometric interpolation setup
2278  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
2279  & geometric_component_number,err,error,*999)
2280 
2281  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2282  & geometric_component_number,err,error,*999)
2283  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
2284  & geometric_component_number,err,error,*999)
2285 
2286  SELECT CASE(equations_set%SOLUTION_METHOD)
2288  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2289  & field_node_based_interpolation,err,error,*999)
2290  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2291  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
2292  !Default the scaling to the geometric field scaling
2293  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
2294  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
2296  CALL flagerror("Not implemented.",err,error,*999)
2298  CALL flagerror("Not implemented.",err,error,*999)
2300  CALL flagerror("Not implemented.",err,error,*999)
2302  CALL flagerror("Not implemented.",err,error,*999)
2304  CALL flagerror("Not implemented.",err,error,*999)
2305  CASE DEFAULT
2306  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2307  & " is invalid."
2308  CALL flagerror(local_error,err,error,*999)
2309  END SELECT
2310  ELSE
2311  !Check the user specified field
2312  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
2313  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
2314  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
2315  CALL field_variable_types_check(equations_set_setup%FIELD, &
2316  & (/field_u_variable_type,field_deludeln_variable_type/),err,error,*999)
2317 
2318  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type,err,error,*999)
2319  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
2320  & err,error,*999)
2321 
2322  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
2323  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
2324 
2325  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
2326  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1,err,error,*999)
2327 
2328  SELECT CASE(equations_set%SOLUTION_METHOD)
2330  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
2331  & field_node_based_interpolation,err,error,*999)
2332  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
2333  & field_node_based_interpolation,err,error,*999)
2335  CALL flagerror("Not implemented.",err,error,*999)
2337  CALL flagerror("Not implemented.",err,error,*999)
2339  CALL flagerror("Not implemented.",err,error,*999)
2341  CALL flagerror("Not implemented.",err,error,*999)
2343  CALL flagerror("Not implemented.",err,error,*999)
2344  CASE DEFAULT
2345  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2346  & " is invalid."
2347  CALL flagerror(local_error,err,error,*999)
2348  END SELECT
2349  ENDIF
2351  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
2352  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
2353  ENDIF
2354  CASE DEFAULT
2355  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2356  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2357  & " is invalid for an extracellular bidomain Poisson equation"
2358  CALL flagerror(local_error,err,error,*999)
2359  END SELECT
2360  !MATERIAL
2362  SELECT CASE(equations_set_setup%ACTION_TYPE)
2363 
2365  equations_materials=>equations_set%MATERIALS
2366  IF(ASSOCIATED(equations_materials)) THEN
2367  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
2368  !Create the auto created materials field
2369  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
2370  & materials_field,err,error,*999)
2371  CALL field_label_set(equations_set%MATERIALS%MATERIALS_FIELD,"Material Field",err,error,*999)
2372  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
2373  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
2374  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
2375  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
2376  & err,error,*999)
2377  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
2378  & geometric_field,err,error,*999)
2379  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
2380  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,(/field_u_variable_type/), &
2381  & err,error,*999)
2382  CALL field_variable_label_set(equations_set%MATERIALS%MATERIALS_FIELD,field_u_variable_type,"conductivity",err, &
2383  & error,*999)
2384  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2385  & field_vector_dimension_type,err,error,*999)
2386  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2387  & field_dp_type,err,error,*999)
2388  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2389  & number_of_dimensions,err,error,*999)
2390  IF(number_of_dimensions==1) THEN
2391  number_of_materials_components=2
2392  ELSEIF(number_of_dimensions==2) THEN
2393  number_of_materials_components=6
2394  ELSEIF(number_of_dimensions==3) THEN
2395  number_of_materials_components=12
2396  ENDIF
2397  !Set the number of materials components
2398  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2399  & number_of_materials_components,err,error,*999)
2400  !Default the materials components to the first geometric component with constant interpolation
2401  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2402  & 1,geometric_component_number,err,error,*999)
2403  DO component_idx=1,number_of_materials_components
2404  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2405  & component_idx,geometric_component_number,err,error,*999)
2406  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2407  & component_idx,field_constant_interpolation,err,error,*999)
2408  ENDDO !components_idx
2409  !Default the field scaling to that of the geometric field
2410  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
2411  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
2412  ELSE
2413  !Check the user specified field
2414  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
2415  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
2416  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
2417  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
2418  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
2419  & err,error,*999)
2420  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
2421  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2422  & number_of_dimensions,err,error,*999)
2423  IF(number_of_dimensions==1) THEN
2424  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,2,err, &
2425  & error,*999)
2426  ELSEIF(number_of_dimensions==2) THEN
2427  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,6,err, &
2428  & error,*999)
2429  ELSEIF(number_of_dimensions==3) THEN
2430  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,12,err, &
2431  & error,*999)
2432  ENDIF
2433  ENDIF
2434  ELSE
2435  CALL flagerror("Equations set materials is not associated.",err,error,*999)
2436  ENDIF
2437 
2439  equations_materials=>equations_set%MATERIALS
2440  IF(ASSOCIATED(equations_materials)) THEN
2441  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
2442  !Finish creating the materials field
2443  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
2444  !Set the default values for the materials field
2445  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2446  & number_of_dimensions,err,error,*999)
2447  IF(number_of_dimensions==1) THEN
2448  number_of_materials_components=2
2449  ELSEIF(number_of_dimensions==2) THEN
2450  number_of_materials_components=6
2451  ELSEIF(number_of_dimensions==3) THEN
2452  number_of_materials_components=12
2453  ENDIF
2454  !Set all material component values to one
2455  DO component_idx=1,number_of_materials_components
2456  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2457  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
2458  ENDDO !component_idx
2459  ENDIF
2460  ELSE
2461  CALL flagerror("Equations set materials is not associated.",err,error,*999)
2462  ENDIF
2463  CASE DEFAULT
2464  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2465  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2466  & " is invalid for an extrcellular bidomain Poisson equation."
2467  CALL flagerror(local_error,err,error,*999)
2468  END SELECT
2469  !SOURCE
2471  SELECT CASE(equations_set%SPECIFICATION(3))
2473  SELECT CASE(equations_set_setup%ACTION_TYPE)
2474 
2475  !Set start action
2477  !Do nothing
2478  IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
2479  !Create the auto created source field
2480  !start field creation with name 'Vm'
2481  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
2482  & equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
2483  !start creation of a new field
2484  CALL field_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_general_type,err,error,*999)
2485  !label the field
2486  CALL field_label_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,"Vm",err,error,*999)
2487  !define new created field to be source
2488  CALL field_dependent_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
2489  & field_independent_type,err,error,*999)
2490  !look for decomposition rule already defined
2491  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
2492  & err,error,*999)
2493  !apply decomposition rule found on new created field
2494  CALL field_mesh_decomposition_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
2495  & geometric_decomposition,err,error,*999)
2496  !point new field to geometric field
2497  CALL field_geometric_field_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,equations_set% &
2498  & geometry%GEOMETRIC_FIELD,err,error,*999)
2499  !set number of variables to 1
2500  source_field_number_of_variables=1
2501  CALL field_number_of_variables_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
2502  & source_field_number_of_variables,err,error,*999)
2503  CALL field_variable_types_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
2504  & (/field_u_variable_type/),err,error,*999)
2505  CALL field_dimension_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
2506  & field_vector_dimension_type,err,error,*999)
2507  CALL field_data_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
2508  & field_dp_type,err,error,*999)
2509  !set the number of components to 1
2510  source_field_number_of_components=1
2511  CALL field_number_of_components_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
2512  & field_u_variable_type,source_field_number_of_components,err,error,*999)
2513  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2514  & 1,geometric_mesh_component,err,error,*999)
2515  !Default to the geometric interpolation setup
2516  CALL field_component_mesh_component_set(equations_set%SOURCE%SOURCE_FIELD, &
2517  & field_u_variable_type,1,geometric_mesh_component,err,error,*999)
2518  SELECT CASE(equations_set%SOLUTION_METHOD)
2519  !Specify fem solution method
2521  CALL field_component_interpolation_set_and_lock(equations_set%SOURCE%SOURCE_FIELD, &
2522  & field_u_variable_type,1,field_node_based_interpolation,err,error,*999)
2523  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
2524  & err,error,*999)
2525  CALL field_scaling_type_set(equations_set%SOURCE%SOURCE_FIELD,geometric_scaling_type, &
2526  & err,error,*999)
2527  !Other solutions not defined yet
2528  CASE DEFAULT
2529  local_error="The solution method of " &
2530  & //trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// " is invalid."
2531  CALL flagerror(local_error,err,error,*999)
2532  END SELECT
2533  ELSE
2534  !Check the user specified field
2535  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
2536  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
2537  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
2538  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
2539  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
2540  & err,error,*999)
2541  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
2542  !number of components
2543  source_field_number_of_components=1
2544  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
2545  & source_field_number_of_components,err,error,*999)
2546  SELECT CASE(equations_set%SOLUTION_METHOD)
2548  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
2549  & field_node_based_interpolation,err,error,*999)
2550  CASE DEFAULT
2551  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
2552  &"*",err,error))//" is invalid."
2553  CALL flagerror(local_error,err,error,*999)
2554  END SELECT
2555  ENDIF
2556 
2557  !Specify finish action
2559  IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
2560  CALL field_create_finish(equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
2561  ENDIF
2562  CASE DEFAULT
2563  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2564  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2565  & " is invalid for an extracellular bidomain Poisson subtype"
2566  CALL flagerror(local_error,err,error,*999)
2567  END SELECT
2568  CASE DEFAULT
2569  local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2570  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2571  & " is invalid for an extracellular bidomain Poisson equation."
2572  CALL flagerror(local_error,err,error,*999)
2573  END SELECT
2574  !ANALYTC
2576  SELECT CASE(equations_set_setup%ACTION_TYPE)
2578  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
2579  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
2580  IF(ASSOCIATED(dependent_field)) THEN
2581  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
2582  IF(ASSOCIATED(geometric_field)) THEN
2583  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
2584  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
2586  !Check that we have an constant source
2587  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
2588  !Check that we are in 2D
2589  IF(number_of_dimensions/=2) THEN
2590  local_error="The number of geometric dimensions of "// &
2591  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2592  & " is invalid. The analytic function type of "// &
2593  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2594  & " requires that there be 2 geometric dimensions."
2595  CALL flagerror(local_error,err,error,*999)
2596  ENDIF
2597  !Create analytic field if required
2598  !Set analtyic function type
2599  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_poisson_equation_two_dim_1
2600  ELSE
2601  local_error="The equations set subtype of "// &
2602  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2603  & " is invalid. The analytic function type of "// &
2604  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2605  & " requires that the equations set subtype be an exponential source Poisson equation."
2606  CALL flagerror(local_error,err,error,*999)
2607  ENDIF
2609  !Check that we have an exponential source
2610  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
2611  !Check that we are in 3D
2612  IF(number_of_dimensions/=3) THEN
2613  local_error="The number of geometric dimensions of "// &
2614  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2615  & " is invalid. The analytic function type of "// &
2616  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2617  & " requires that there be 3 geometric dimensions."
2618  CALL flagerror(local_error,err,error,*999)
2619  ENDIF
2620  !Create analytic field if required
2621  !Set analytic function type
2622  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_poisson_equation_three_dim_1
2623  ELSE
2624  local_error="The equations set subtype of "// &
2625  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2626  & " is invalid. The analytic function type of "// &
2627  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2628  & " requires that the equations set subtype be an exponential source Poisson equation."
2629  CALL flagerror(local_error,err,error,*999)
2630  ENDIF
2631  CASE DEFAULT
2632  local_error="The specified analytic function type of "// &
2633  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2634  & " is invalid for a linear source Poisson equation."
2635  CALL flagerror(local_error,err,error,*999)
2636  END SELECT
2637  ELSE
2638  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
2639  ENDIF
2640  ELSE
2641  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
2642  ENDIF
2643  ELSE
2644  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2645  ENDIF
2647  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
2648  analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
2649  IF(ASSOCIATED(analytic_field)) THEN
2650  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
2651  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
2652  ENDIF
2653  ENDIF
2654  ELSE
2655  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
2656  ENDIF
2657  CASE DEFAULT
2658  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2659  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2660  & " is invalid for a extracellular bidomain source Poisson equation."
2661  CALL flagerror(local_error,err,error,*999)
2662  END SELECT
2663  !EQUATIONS
2665  SELECT CASE(equations_set_setup%ACTION_TYPE)
2666  !start action
2668  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
2669  !Create the equations
2670  CALL equations_create_start(equations_set,equations,err,error,*999)
2671  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
2672  CALL equations_time_dependence_type_set(equations,equations_static,err,error,*999)
2673  ELSE
2674  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2675  ENDIF
2676  !finish action
2678  SELECT CASE(equations_set%SOLUTION_METHOD)
2680  !Finish the creation of the equations
2681  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
2682  CALL equations_create_finish(equations,err,error,*999)
2683  !Create the equations mapping.
2684  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
2685  IF(equations_set%SPECIFICATION(3)==equations_set_extracellular_bidomain_poisson_subtype) THEN
2686  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
2687  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,[field_u_variable_type],err,error,*999)
2688  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
2689  CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
2690  ELSE
2691  local_error="The third equations set specification of "// &
2692  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))//" is invalid."
2693  CALL flagerror(local_error,err,error,*999)
2694  ENDIF
2695  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
2696  !Create the equations matrices
2697  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
2698  SELECT CASE(equations%SPARSITY_TYPE)
2700  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
2701  & (/matrix_block_storage_type/),err,error,*999)
2703  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
2704  & (/matrix_compressed_row_storage_type/),err,error,*999)
2705  CALL equationsmatrices_linearstructuretypeset(equations_matrices, &
2706  & (/equations_matrix_fem_structure/),err,error,*999)
2707  CASE DEFAULT
2708  local_error="The equations matrices sparsity type of "// &
2709  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
2710  CALL flagerror(local_error,err,error,*999)
2711  END SELECT
2712  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
2714  CALL flagerror("Not implemented.",err,error,*999)
2716  CALL flagerror("Not implemented.",err,error,*999)
2718  CALL flagerror("Not implemented.",err,error,*999)
2720  CALL flagerror("Not implemented.",err,error,*999)
2722  CALL flagerror("Not implemented.",err,error,*999)
2723  CASE DEFAULT
2724  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2725  & " is invalid."
2726  CALL flagerror(local_error,err,error,*999)
2727  END SELECT
2728  CASE DEFAULT
2729  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2730  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2731  & " is invalid for a linear source Poisson equation."
2732  CALL flagerror(local_error,err,error,*999)
2733  END SELECT
2734  CASE DEFAULT
2735  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2736  & " is invalid for a extracellular bidomain source Poisson equation."
2737  CALL flagerror(local_error,err,error,*999)
2738  END SELECT
2739  ELSE
2740  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2741  & " is not a extracellular bidomain source Poisson equation subtype."
2742  CALL flagerror(local_error,err,error,*999)
2743  ENDIF
2744  ELSE
2745  CALL flagerror("Equations set is not associated.",err,error,*999)
2746  ENDIF
2747 
2748  exits("Poisson_EquationsSetExtracellularBidomainSetup")
2749  RETURN
2750 999 errors("Poisson_EquationsSetExtracellularBidomainSetup",err,error)
2751  exits("Poisson_EquationsSetExtracellularBidomainSetup")
2752  RETURN 1
2753 
2755 
2756  !
2757  !================================================================================================================================
2758  !
2759 
2761  SUBROUTINE poisson_equationssetnonlinearsourcesetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
2763  !Argument variables
2764  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2765  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
2766  INTEGER(INTG), INTENT(OUT) :: ERR
2767  TYPE(varying_string), INTENT(OUT) :: ERROR
2768  !Local Variables
2769  INTEGER(INTG) :: component_idx,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE,NUMBER_OF_DIMENSIONS, &
2770  & NUMBER_OF_MATERIALS_COMPONENTS
2771  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
2772  TYPE(equations_type), POINTER :: EQUATIONS
2773  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
2774  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
2775  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
2776  TYPE(field_type), POINTER :: analytic_field,DEPENDENT_FIELD,geometric_field
2777  TYPE(varying_string) :: LOCAL_ERROR
2778 
2779  enters("POISSON_EQUATION_EQUATION_SET_NONLINEAR_SOURCE_SETUP",err,error,*999)
2780 
2781  NULLIFY(equations)
2782  NULLIFY(equations_mapping)
2783  NULLIFY(equations_matrices)
2784  NULLIFY(geometric_decomposition)
2785 
2786  IF(ASSOCIATED(equations_set)) THEN
2787  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
2788  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
2789  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
2790  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
2791  & err,error,*999)
2792  END IF
2793  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_poisson_subtype.OR. &
2794  & equations_set%SPECIFICATION(3)==equations_set_exponential_source_poisson_subtype) THEN
2795  SELECT CASE(equations_set_setup%SETUP_TYPE)
2797  SELECT CASE(equations_set_setup%ACTION_TYPE)
2801 !!TODO: Check valid setup
2802  CASE DEFAULT
2803  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2804  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2805  & " is invalid for a nonlinear source Poisson equation."
2806  CALL flagerror(local_error,err,error,*999)
2807  END SELECT
2809  !Do nothing???
2811  SELECT CASE(equations_set_setup%ACTION_TYPE)
2813  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
2814  !Create the auto created dependent field
2815  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
2816  & dependent_field,err,error,*999)
2817  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
2818  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
2819  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
2820  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
2821  & err,error,*999)
2822  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
2823  & geometric_field,err,error,*999)
2824  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
2825  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,(/field_u_variable_type, &
2826  & field_deludeln_variable_type/),err,error,*999)
2827  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2828  & field_scalar_dimension_type,err,error,*999)
2829  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
2830  & field_scalar_dimension_type,err,error,*999)
2831  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2832  & field_dp_type,err,error,*999)
2833  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
2834  & field_dp_type,err,error,*999)
2835  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2836  & err,error,*999)
2837  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
2838  & err,error,*999)
2839  !Default to the geometric interpolation setup
2840  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
2841  & geometric_component_number,err,error,*999)
2842  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2843  & geometric_component_number,err,error,*999)
2844  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
2845  & geometric_component_number,err,error,*999)
2846  SELECT CASE(equations_set%SOLUTION_METHOD)
2848  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2849  & field_node_based_interpolation,err,error,*999)
2850  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2851  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
2852  !Default the scaling to the geometric field scaling
2853  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
2854  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
2856  CALL flagerror("Not implemented.",err,error,*999)
2858  CALL flagerror("Not implemented.",err,error,*999)
2860  CALL flagerror("Not implemented.",err,error,*999)
2862  CALL flagerror("Not implemented.",err,error,*999)
2864  CALL flagerror("Not implemented.",err,error,*999)
2865  CASE DEFAULT
2866  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2867  & " is invalid."
2868  CALL flagerror(local_error,err,error,*999)
2869  END SELECT
2870  ELSE
2871  !Check the user specified field
2872  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
2873  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
2874  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
2875  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type,field_deludeln_variable_type/), &
2876  & err,error,*999)
2877  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type,err,error,*999)
2878  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
2879  & err,error,*999)
2880  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
2881  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
2882  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
2883  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1,err,error,*999)
2884  SELECT CASE(equations_set%SOLUTION_METHOD)
2886  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
2887  & field_node_based_interpolation,err,error,*999)
2888  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
2889  & field_node_based_interpolation,err,error,*999)
2891  CALL flagerror("Not implemented.",err,error,*999)
2893  CALL flagerror("Not implemented.",err,error,*999)
2895  CALL flagerror("Not implemented.",err,error,*999)
2897  CALL flagerror("Not implemented.",err,error,*999)
2899  CALL flagerror("Not implemented.",err,error,*999)
2900  CASE DEFAULT
2901  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2902  & " is invalid."
2903  CALL flagerror(local_error,err,error,*999)
2904  END SELECT
2905  ENDIF
2907  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
2908  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
2909  ENDIF
2910  CASE DEFAULT
2911  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2912  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2913  & " is invalid for a nonlinear source Poisson equation"
2914  CALL flagerror(local_error,err,error,*999)
2915  END SELECT
2917  SELECT CASE(equations_set_setup%ACTION_TYPE)
2919  equations_materials=>equations_set%MATERIALS
2920  IF(ASSOCIATED(equations_materials)) THEN
2921  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
2922  !Create the auto created materials field
2923  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
2924  & materials_field,err,error,*999)
2925  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
2926  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
2927  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
2928  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
2929  & err,error,*999)
2930  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
2931  & geometric_field,err,error,*999)
2932  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
2933  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2934  & field_vector_dimension_type,err,error,*999)
2935  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2936  & field_dp_type,err,error,*999)
2937  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2938  & number_of_dimensions,err,error,*999)
2939  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_poisson_subtype) THEN
2940  !Quadratic source. Materials field components are 1 for each dimension and 3 for the quadratic source
2941  !i.e., k and a, b and c in div(k.grad(u(x)))=a(x)+b(x)u(x)+c(x)u^2(x)
2942  number_of_materials_components=number_of_dimensions+3
2943  ELSE
2944  !Exponential source. Matierals field components are 1 for each dimension and 3 for the exponential source
2945  !i.e., k, a, b and c in div(k.grad(u(x)))=a(x)+b(x)e^[c(x)u(x)]
2946  number_of_materials_components=number_of_dimensions+3
2947  ENDIF
2948  !Set the number of materials components
2949  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2950  & number_of_materials_components,err,error,*999)
2951  !Default the k materials components to the geometric interpolation setup with constant interpolation
2952  DO component_idx=1,number_of_dimensions
2953  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2954  & component_idx,geometric_component_number,err,error,*999)
2955  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2956  & component_idx,geometric_component_number,err,error,*999)
2957  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2958  & component_idx,field_constant_interpolation,err,error,*999)
2959  ENDDO !component_idx
2960  !Default the source materials components to the first component geometric interpolation with constant interpolation
2961  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2962  & 1,geometric_component_number,err,error,*999)
2963  DO component_idx=number_of_dimensions+1,number_of_materials_components
2964  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2965  & component_idx,geometric_component_number,err,error,*999)
2966  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2967  & component_idx,field_constant_interpolation,err,error,*999)
2968  ENDDO !components_idx
2969  !Default the field scaling to that of the geometric field
2970  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
2971  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
2972  ELSE
2973  !Check the user specified field
2974  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
2975  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
2976  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
2977  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
2978  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
2979  & err,error,*999)
2980  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
2981  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2982  & number_of_dimensions,err,error,*999)
2983  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_poisson_subtype) THEN
2984  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+3, &
2985  & err,error,*999)
2986  ELSE
2987  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+3, &
2988  & err,error,*999)
2989  ENDIF
2990  ENDIF
2991  ELSE
2992  CALL flagerror("Equations set materials is not associated.",err,error,*999)
2993  ENDIF
2995  equations_materials=>equations_set%MATERIALS
2996  IF(ASSOCIATED(equations_materials)) THEN
2997  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
2998  !Finish creating the materials field
2999  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
3000  !Set the default values for the materials field
3001  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
3002  & number_of_dimensions,err,error,*999)
3003  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_poisson_subtype) THEN
3004  !Quadratic source. Materials field components are 1 for each dimension and 3 for the quadratic source
3005  !i.e., k and a, b and c in div(k.grad(u(x)))=a(x)+b(x)u(x)+c(x)u^2(x)
3006  number_of_materials_components=number_of_dimensions+3
3007  ELSE
3008  !Exponential source. Matierals field components are 1 for each dimension and 3 for the exponential source
3009  !i.e., k, a, b and c in div(k.grad(u(x)))=a(x)+b(x)e^[c(x)u(x)]
3010  number_of_materials_components=number_of_dimensions+3
3011  ENDIF
3012  !First set the k values to 1.0
3013  DO component_idx=1,number_of_dimensions
3014  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
3015  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
3016  ENDDO !component_idx
3017  !Now set the source values to 1.0
3018  DO component_idx=number_of_dimensions+1,number_of_materials_components
3019  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
3020  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
3021  ENDDO !component_idx
3022  ENDIF
3023  ELSE
3024  CALL flagerror("Equations set materials is not associated.",err,error,*999)
3025  ENDIF
3026  CASE DEFAULT
3027  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
3028  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
3029  & " is invalid for a linear source Poisson equation."
3030  CALL flagerror(local_error,err,error,*999)
3031  END SELECT
3033  SELECT CASE(equations_set_setup%ACTION_TYPE)
3035  !Do nothing
3037  !Do nothing
3038  !? Maybe set finished flag????
3039  CASE DEFAULT
3040  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
3041  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
3042  & " is invalid for a nonlinear source Poisson equation."
3043  CALL flagerror(local_error,err,error,*999)
3044  END SELECT
3046  SELECT CASE(equations_set_setup%ACTION_TYPE)
3048  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
3049  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
3050  IF(ASSOCIATED(dependent_field)) THEN
3051  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
3052  IF(ASSOCIATED(geometric_field)) THEN
3053  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
3054  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
3056  !Check that we have an exponential source
3057  IF(equations_set%SPECIFICATION(3)==equations_set_exponential_source_poisson_subtype) THEN
3058  !Check that we are in 2D
3059  IF(number_of_dimensions/=2) THEN
3060  local_error="The number of geometric dimensions of "// &
3061  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
3062  & " is invalid. The analytic function type of "// &
3063  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
3064  & " requires that there be 2 geometric dimensions."
3065  CALL flagerror(local_error,err,error,*999)
3066  ENDIF
3067  !Create analytic field if required
3068  !Set analytic function type
3069  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_poisson_equation_two_dim_1
3070  ELSE
3071  local_error="The equations set subtype of "// &
3072  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
3073  & " is invalid. The analytic function type of "// &
3074  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
3075  & " requires that the equations set subtype be an exponential source Poisson equation."
3076  CALL flagerror(local_error,err,error,*999)
3077  ENDIF
3079  !Check that we have an exponential source
3080  IF(equations_set%SPECIFICATION(3)==equations_set_exponential_source_poisson_subtype) THEN
3081  !Check that we are in 3D
3082  IF(number_of_dimensions/=3) THEN
3083  local_error="The number of geometric dimensions of "// &
3084  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
3085  & " is invalid. The analytic function type of "// &
3086  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
3087  & " requires that there be 3 geometric dimensions."
3088  CALL flagerror(local_error,err,error,*999)
3089  ENDIF
3090  !Create analytic field if required
3091  !Set analytic function type
3092  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_poisson_equation_three_dim_1
3093  ELSE
3094  local_error="The equations set subtype of "// &
3095  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
3096  & " is invalid. The analytic function type of "// &
3097  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
3098  & " requires that the equations set subtype be an exponential source Poisson equation."
3099  CALL flagerror(local_error,err,error,*999)
3100  ENDIF
3101  CASE DEFAULT
3102  local_error="The specified analytic function type of "// &
3103  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
3104  & " is invalid for a nonlinear source Poisson equation."
3105  CALL flagerror(local_error,err,error,*999)
3106  END SELECT
3107  ELSE
3108  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
3109  ENDIF
3110  ELSE
3111  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
3112  ENDIF
3113  ELSE
3114  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
3115  ENDIF
3117  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
3118  analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
3119  IF(ASSOCIATED(analytic_field)) THEN
3120  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
3121  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
3122  ENDIF
3123  ENDIF
3124  ELSE
3125  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
3126  ENDIF
3127  CASE DEFAULT
3128  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
3129  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
3130  & " is invalid for a nonlinear source Poisson equation."
3131  CALL flagerror(local_error,err,error,*999)
3132  END SELECT
3134  SELECT CASE(equations_set_setup%ACTION_TYPE)
3136  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
3137  !Start the creation of the equations
3138  CALL equations_create_start(equations_set,equations,err,error,*999)
3139  CALL equations_linearity_type_set(equations,equations_nonlinear,err,error,*999)
3140  CALL equations_time_dependence_type_set(equations,equations_static,err,error,*999)
3141  ELSE
3142  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
3143  ENDIF
3145  SELECT CASE(equations_set%SOLUTION_METHOD)
3147  !Finish the creation of the equations
3148  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
3149  CALL equations_create_finish(equations,err,error,*999)
3150  !Create the equations mapping.
3151  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
3152  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
3153  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,(/field_u_variable_type/), &
3154  & err,error,*999)
3155  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
3156  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
3157  !Create the equations matrices
3158  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
3159  ! Use the analytic Jacobian calculation
3161  & err,error,*999)
3162  SELECT CASE(equations%SPARSITY_TYPE)
3165  & err,error,*999)
3167  & err,error,*999)
3170  & err,error,*999)
3172  & err,error,*999)
3174  & err,error,*999)
3176  & err,error,*999)
3177  CASE DEFAULT
3178  local_error="The equations matrices sparsity type of "// &
3179  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
3180  CALL flagerror(local_error,err,error,*999)
3181  END SELECT
3182  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
3184  CALL flagerror("Not implemented.",err,error,*999)
3186  CALL flagerror("Not implemented.",err,error,*999)
3188  CALL flagerror("Not implemented.",err,error,*999)
3190  CALL flagerror("Not implemented.",err,error,*999)
3192  CALL flagerror("Not implemented.",err,error,*999)
3193  CASE DEFAULT
3194  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
3195  & " is invalid."
3196  CALL flagerror(local_error,err,error,*999)
3197  END SELECT
3198  CASE DEFAULT
3199  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
3200  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
3201  & " is invalid for a nonlinear source Poisson equation."
3202  CALL flagerror(local_error,err,error,*999)
3203  END SELECT
3204  CASE DEFAULT
3205  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
3206  & " is invalid for a nonlinear source Poisson equation."
3207  CALL flagerror(local_error,err,error,*999)
3208  END SELECT
3209  ELSE
3210  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
3211  & " is not a nonlinear source Poisson equation subtype."
3212  CALL flagerror(local_error,err,error,*999)
3213  ENDIF
3214  ELSE
3215  CALL flagerror("Equations set is not associated.",err,error,*999)
3216  ENDIF
3217 
3218  exits("Poisson_EquationsSetNonlinearSourceSetup")
3219  RETURN
3220 999 errors("Poisson_EquationsSetNonlinearSourceSetup",err,error)
3221  exits("Poisson_EquationsSetNonlinearSourceSetup")
3222  RETURN 1
3223 
3225 
3226  !
3227  !================================================================================================================================
3228  !
3229 
3231  SUBROUTINE poisson_equation_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
3233  !Argument variables
3234  TYPE(problem_type), POINTER :: PROBLEM
3235  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
3236  INTEGER(INTG), INTENT(OUT) :: ERR
3237  TYPE(varying_string), INTENT(OUT) :: ERROR
3238  !Local Variables
3239  TYPE(varying_string) :: LOCAL_ERROR
3240 
3241  enters("POISSON_EQUATION_PROBLEM_SETUP",err,error,*999)
3242 
3243  IF(ASSOCIATED(problem)) THEN
3244  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
3245  CALL flagerror("Problem specification is not allocated.",err,error,*999)
3246  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
3247  CALL flagerror("Problem specification must have three entries for a Poisson problem.",err,error,*999)
3248  END IF
3249  SELECT CASE(problem%SPECIFICATION(3))
3251  CALL poisson_problemextracellularbidomainsetup(problem,problem_setup,err,error,*999)
3253  CALL poisson_problemlinearsourcesetup(problem,problem_setup,err,error,*999)
3256  CALL poisson_problempressurepoissonsetup(problem,problem_setup,err,error,*999)
3258  CALL poisson_problemnonlinearsourcesetup(problem,problem_setup,err,error,*999)
3259  CASE DEFAULT
3260  local_error="Problem subtype "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
3261  & " is not valid for a Poisson equation type of a classical field problem class."
3262  CALL flagerror(local_error,err,error,*999)
3263  END SELECT
3264  ELSE
3265  CALL flagerror("Problem is not associated.",err,error,*999)
3266  ENDIF
3267 
3268  exits("POISSON_EQUATION_PROBLEM_SETUP")
3269  RETURN
3270 999 errorsexits("POISSON_EQUATION_PROBLEM_SETUP",err,error)
3271  RETURN 1
3272  END SUBROUTINE poisson_equation_problem_setup
3273 
3274  !
3275  !================================================================================================================================
3276  !
3277 
3279  SUBROUTINE poisson_equation_finite_element_calculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
3281  !Argument variables
3282  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
3283  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
3284  INTEGER(INTG), INTENT(OUT) :: ERR
3285  TYPE(varying_string), INTENT(OUT) :: ERROR
3286  !Local Variables
3287  INTEGER(INTG) FIELD_VAR_TYPE,ng,mh,mhs,mi,ms,nh,nhs,ni,ns,I,J,K,L,H,element_node_identity !,k1sum,k2sum
3288  REAL(DP) :: RWG,SUM,PGMSI(3),PGNSI(3),SUM2,DXI_DX(3,3),DELTA_T,DXI_DX_DX(3,3),PHINS
3289  REAL(DP) :: CONDUCTIVITY_SUM_MATERIAL(3,3),CONDUCTIVITY_SUM(3,3),CONDUCTIVITY_SUM_TEMP(3,3),dNudX(3,3),dXdNu(3,3), &
3290  & DNUDXI(3,3),DXIDNU(3,3)
3291  REAL(DP) :: CONDUCTIVITY_I_MATERIAL(3,3),CONDUCTIVITY_I(3,3),CONDUCTIVITY_I_TEMP(3,3)
3292  REAL(DP), ALLOCATABLE :: MATRIX_K_I(:,:)
3293  REAL(DP), ALLOCATABLE :: Vm(:)
3294  REAL(DP) :: U_VALUE(3),U_DERIV(3,3),RHO_PARAM,MU_PARAM,U_OLD(3),U_SECOND(3,3,3),X(3),B(3),P_DERIV(3),W_VALUE(3)
3295  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS,SOURCE_BASIS,INDEPENDENT_BASIS
3296  TYPE(equations_type), POINTER :: EQUATIONS
3297  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
3298  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
3299  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
3300  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
3301  TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
3302  TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
3303  TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
3304  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,SOURCE_FIELD,MATERIALS_FIELD,INDEPENDENT_FIELD,FIBRE_FIELD
3305  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,SOURCE_FIELD_VARIABLE
3306  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
3307  TYPE(field_interpolated_point_type), POINTER :: GEOMETRIC_INTERPOLATED_POINT,FIBRE_INTERPOLATED_POINT
3308  TYPE(field_interpolated_point_metrics_type), POINTER :: GEOMETRIC_INTERP_POINT_METRICS
3309  TYPE(varying_string) :: LOCAL_ERROR
3310 
3311  INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,NUMBER_OF_XI,NUMBER_OF_ELEMENT_PARAMETERS,node_idx,dof_idx
3312 
3313  LOGICAL :: INSIDE,BETWEEN
3314  REAL(DP), POINTER :: INPUT_LABEL(:)
3315  REAL(DP) :: DIFF_COEFF1,DIFF_COEFF2
3316 
3317 #ifdef TAUPROF
3318  CHARACTER(26) :: CVAR
3319  INTEGER :: GAUSS_POINT_LOOP_PHASE(2) = [ 0, 0 ]
3320  SAVE gauss_point_loop_phase
3321 #endif
3322 
3323  NULLIFY(source_field_variable)
3324 
3325  enters("POISSON_EQUATION_FINITE_ELEMENT_CALCULATE",err,error,*999)
3326 
3327  IF(ASSOCIATED(equations_set)) THEN
3328  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
3329  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
3330  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
3331  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
3332  & err,error,*999)
3333  END IF
3334  equations=>equations_set%EQUATIONS
3335  IF(ASSOCIATED(equations)) THEN
3336  SELECT CASE(equations_set%SPECIFICATION(3))
3338  !Store all these in equations matrices/somewhere else?????
3339  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
3340  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
3341  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
3342  equations_matrices=>equations%EQUATIONS_MATRICES
3343  linear_matrices=>equations_matrices%LINEAR_MATRICES
3344  equations_matrix=>linear_matrices%MATRICES(1)%PTR
3345  rhs_vector=>equations_matrices%RHS_VECTOR
3346  source_vector=>equations_matrices%SOURCE_VECTOR
3347  equations_mapping=>equations%EQUATIONS_MAPPING
3348  linear_mapping=>equations_mapping%LINEAR_MAPPING
3349  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
3350  field_var_type=field_variable%VARIABLE_TYPE
3351  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3352  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3353  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3354  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3355  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
3356  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3357  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3358  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3359  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3360  !Loop over gauss points
3361  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
3362  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3363  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
3364  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3365  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
3366  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
3367  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
3368  !Calculate RWG.
3369 !!TODO: Think about symmetric problems.
3370  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
3371  & quadrature_scheme%GAUSS_WEIGHTS(ng)
3372  !Loop over field components
3373  mhs=0
3374  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3375  !Loop over element rows
3376 !!TODO: CHANGE ELEMENT CALCULATE TO WORK OF ns ???
3377  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
3378  mhs=mhs+1
3379  nhs=0
3380  IF(equations_matrix%UPDATE_MATRIX) THEN
3381  !Loop over element columns
3382  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
3383  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
3384  nhs=nhs+1
3385  DO ni=1,dependent_basis%NUMBER_OF_XI
3386  pgmsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
3387  pgnsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
3388  ENDDO !ni
3389  sum=0.0_dp
3390  DO mi=1,dependent_basis%NUMBER_OF_XI
3391  DO ni=1,dependent_basis%NUMBER_OF_XI
3392  sum=sum+pgmsi(mi)*pgnsi(ni)*equations%INTERPOLATION% &
3393  & geometric_interp_point_metrics(field_u_variable_type)%PTR%GU(mi,ni)
3394  ENDDO !ni
3395  ENDDO !mi
3396  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg
3397  ENDDO !ns
3398  ENDDO !nh
3399  ENDIF
3400  IF(source_vector%UPDATE_VECTOR) THEN
3401  !Use materials field value
3402  sum=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
3403  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+ & !This needs to changed to the source vector
3404  & sum*quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*rwg
3405  ENDIF
3406  ENDDO !ms
3407  ENDDO !mh
3408  ENDDO !ng
3409 
3411 
3412  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
3413  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
3414  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
3415  fibre_field=>equations%INTERPOLATION%FIBRE_FIELD
3416  source_field=>equations%INTERPOLATION%SOURCE_FIELD
3417 
3418  equations_matrices=>equations%EQUATIONS_MATRICES
3419  linear_matrices=>equations_matrices%LINEAR_MATRICES
3420  equations_matrix=>linear_matrices%MATRICES(1)%PTR
3421 
3422  rhs_vector=>equations_matrices%RHS_VECTOR
3423  source_vector=>equations_matrices%SOURCE_VECTOR
3424  source_vector%UPDATE_VECTOR=.true. !TODO -- maybe done somewhere else (greped at equations_matrices_routines.f90)?, check when time dependency is included
3425 
3426  equations_mapping=>equations%EQUATIONS_MAPPING
3427  linear_mapping=>equations_mapping%LINEAR_MAPPING
3428  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
3429  field_var_type=field_variable%VARIABLE_TYPE
3430 
3431  CALL field_variable_get(source_field,field_u_variable_type,source_field_variable,err,error,*999)
3432 
3433  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3434  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3435  source_basis=>source_field%DECOMPOSITION%DOMAIN(source_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3436  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3437  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3438  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3439 
3440  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
3441  number_of_element_parameters=dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
3442 
3443  number_of_dimensions=equations_set%REGION%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
3444  number_of_xi=dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY% &
3445  & elements%ELEMENTS(element_number)%BASIS%NUMBER_OF_XI
3446 
3447  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3448  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3449  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3450  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3451  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3452  & fibre_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3453  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3454  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
3455  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3456  & source_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3457 
3458  geometric_interpolated_point=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
3459  fibre_interpolated_point=>equations%INTERPOLATION%FIBRE_INTERP_POINT(field_u_variable_type)%PTR
3460  geometric_interp_point_metrics=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR
3461 
3462  ALLOCATE(matrix_k_i(number_of_element_parameters,number_of_element_parameters),stat=err)
3463  IF(err/=0) CALL flagerror("Could not allocate matrix K_i.",err,error,*999)
3464  matrix_k_i=0.0_dp
3465  ALLOCATE(vm(number_of_element_parameters),stat=err)
3466  IF(err/=0) CALL flagerror("Could not allocate vector V_m.",err,error,*999)
3467  vm=0.0_dp
3468 
3469  !Loop over gauss points
3470  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
3471 #ifdef TAUPROF
3472  WRITE (cvar,'(a17,i2)') 'Gauss Point Loop ',ng
3473  CALL tau_phase_create_dynamic(gauss_point_loop_phase,cvar)
3474  CALL tau_phase_start(gauss_point_loop_phase)
3475 #endif
3476  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3477  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
3478  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3479  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
3480  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3481  & fibre_interp_point(field_u_variable_type)%PTR,err,error,*999)
3482  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3483  & source_interp_point(field_u_variable_type)%PTR,err,error,*999)
3484 
3485  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,geometric_interp_point_metrics, &
3486  & err,error,*999)
3487 
3488  !Calculate RWG.
3489  rwg=geometric_interp_point_metrics%JACOBIAN*quadrature_scheme%GAUSS_WEIGHTS(ng)
3490  !conductivities in material coordinates
3491  conductivity_sum_material=0.0_dp !intra- plus extracellular
3492  conductivity_i_material=0.0_dp !intracellular
3493  IF(number_of_dimensions==2) THEN
3494  !interpolated values
3495  conductivity_i_material(1,1)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,1)
3496  conductivity_i_material(2,2)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,1)
3497  conductivity_i_material(1,2)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,1)
3498  conductivity_i_material(2,1)=conductivity_i_material(1,2)
3499 
3500  conductivity_sum_material(1,1)=conductivity_i_material(1,1)+ &
3501  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(4,1)
3502  conductivity_sum_material(2,2)=conductivity_i_material(2,2)+ &
3503  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(5,1)
3504  conductivity_sum_material(1,2)=conductivity_i_material(1,2)+ &
3505  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(6,1)
3506  conductivity_sum_material(2,1)=conductivity_sum_material(1,2)
3507  ELSE
3508  !interpolated values
3509  conductivity_i_material(1,1)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,1)
3510  conductivity_i_material(2,2)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,1)
3511  conductivity_i_material(3,3)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,1)
3512  conductivity_i_material(1,2)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(4,1)
3513  conductivity_i_material(2,1)=conductivity_i_material(1,2)
3514  conductivity_i_material(2,3)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(5,1)
3515  conductivity_i_material(3,2)=conductivity_i_material(2,3)
3516  conductivity_i_material(1,3)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(6,1)
3517  conductivity_i_material(3,1)=conductivity_i_material(1,3)
3518 
3519  conductivity_sum_material(1,1)=conductivity_i_material(1,1)+ &
3520  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(7,1)
3521  conductivity_sum_material(2,2)=conductivity_i_material(2,2)+ &
3522  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(8,1)
3523  conductivity_sum_material(3,3)=conductivity_i_material(3,3)+ &
3524  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(9,1)
3525  conductivity_sum_material(1,2)=conductivity_i_material(1,2)+ &
3526  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(10,1)
3527  conductivity_sum_material(2,1)=conductivity_sum_material(1,2)
3528  conductivity_sum_material(2,3)=conductivity_i_material(2,3)+ &
3529  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(11,1)
3530  conductivity_sum_material(3,2)=conductivity_sum_material(2,3)
3531  conductivity_sum_material(1,3)=conductivity_i_material(1,3)+ &
3532  & equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(12,1)
3533  conductivity_sum_material(3,1)=conductivity_sum_material(1,3)
3534  ENDIF
3535 
3536  !rotate the conductivities from material coordinates (nu) to xi-space to get the effective conductivity
3537  CALL coordinates_materialsystemcalculate(geometric_interp_point_metrics,fibre_interpolated_point, &
3538  & dnudx,dxdnu,dnudxi,dxidnu,err,error,*999)
3539 
3540  CALL matrix_product(dnudxi,conductivity_sum_material,conductivity_sum_temp,err,error,*999)
3541  CALL matrix_product(conductivity_sum_temp,dxidnu,conductivity_sum,err,error,*999)
3542 
3543  CALL matrix_product(dnudxi,conductivity_i_material,conductivity_i_temp,err,error,*999)
3544  CALL matrix_product(conductivity_i_temp,dxidnu,conductivity_i,err,error,*999)
3545 
3546  !calculate the element stiffness matrix K_{i+e}
3547  IF(equations_matrix%UPDATE_MATRIX) THEN
3548  !Loop over field components
3549  mhs=0
3550  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3551  !Loop over element rows
3552  DO ms=1,number_of_element_parameters
3553  mhs=mhs+1
3554  !Loop over field components
3555  nhs=0
3556  !Loop over element columns
3557  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
3558  DO ns=1,number_of_element_parameters
3559  nhs=nhs+1
3560 
3561  DO ni=1,dependent_basis%NUMBER_OF_XI
3562  pgmsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
3563  pgnsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
3564  ENDDO !ni
3565 
3566  sum=0.0_dp
3567  DO i=1,dependent_basis%NUMBER_OF_XI
3568  DO k=1,dependent_basis%NUMBER_OF_XI
3569  DO h=1,dependent_basis%NUMBER_OF_XI
3570  sum=sum+conductivity_sum(i,k)*pgnsi(k)*pgmsi(h)*geometric_interp_point_metrics%GU(i,h)
3571  ENDDO !H
3572  ENDDO !K
3573  ENDDO !I
3574 
3575  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg
3576 
3577  ENDDO !ns
3578  ENDDO !nh
3579  ENDDO !ms
3580  ENDDO !mh
3581  ENDIF
3582 
3583  !calculate the matrix K_i for the source vector K_i * Vm
3584  IF(source_vector%UPDATE_VECTOR) THEN
3585  !Loop over field components
3586  mhs=0
3587  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3588  !Loop over element rows
3589  DO ms=1,number_of_element_parameters
3590  mhs=mhs+1
3591  !Loop over field components
3592  nhs=0
3593  !Loop over element columns
3594  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
3595  DO ns=1,number_of_element_parameters
3596  nhs=nhs+1
3597 
3598  DO ni=1,dependent_basis%NUMBER_OF_XI
3599  pgmsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
3600  pgnsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
3601  ENDDO !ni
3602 
3603  sum=0.0_dp
3604  DO i=1,dependent_basis%NUMBER_OF_XI
3605  DO k=1,dependent_basis%NUMBER_OF_XI
3606  DO h=1,dependent_basis%NUMBER_OF_XI
3607  sum=sum+conductivity_i(i,k)*pgnsi(k)*pgmsi(h)*geometric_interp_point_metrics%GU(i,h)
3608  ENDDO !H
3609  ENDDO !K
3610  ENDDO !I
3611 
3612  matrix_k_i(mhs,nhs)=matrix_k_i(mhs,nhs)+sum*rwg
3613 
3614  ENDDO !ns
3615  ENDDO !nh
3616  ENDDO !ms
3617  ENDDO !mh
3618  ENDIF
3619 
3620  ! initialise element RHS_VECTORs with 0
3621  IF(rhs_vector%UPDATE_VECTOR) THEN
3622  !Loop over field components
3623  mhs=0
3624  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3625  !Loop over vector rows
3626  DO ms=1,number_of_element_parameters
3627  mhs=mhs+1
3628  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
3629  ENDDO !ms
3630  ENDDO !mh
3631  ENDIF
3632 
3633 
3634 #ifdef TAUPROF
3635  CALL tau_phase_stop(gauss_point_loop_phase)
3636 #endif
3637  ENDDO !ng
3638 
3639  !now do a matrix-vector of the K_i matrix with the source field Vm and store in source_vector
3640  IF(source_vector%UPDATE_VECTOR) THEN
3641 
3642  !get correct Vm-vector entries (rows nhs)
3643  DO nhs=1,number_of_element_parameters
3644  !get correct global node index for this element
3645  node_idx=source_vector%ELEMENT_VECTOR%ROW_DOFS(nhs)
3646  !get corresponding dof_idx (in this case node_idx=dof_idx as just 1 component, but leave for generality)
3647  dof_idx=source_field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
3648  & derivatives(1)%VERSIONS(1)
3649  !get the correct value from global source field
3650  CALL field_parameter_set_get_local_dof(source_field,field_u_variable_type,field_values_set_type, &
3651  & dof_idx,vm(nhs),err,error,*999)
3652  ENDDO
3653 
3654  !calculate the matrix-vector-product
3655  mhs=0
3656  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3657  !Loop over element matrix / source-vector rows
3658  DO ms=1,number_of_element_parameters
3659  mhs=mhs+1
3660  nhs=0
3661  !Loop over element matrix columns / Vm-vector entries (summation)
3662  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
3663  DO ns=1,number_of_element_parameters
3664  nhs=nhs+1
3665  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+matrix_k_i(mhs,nhs)*vm(nhs)
3666  ENDDO !ns
3667  ENDDO !nh
3668  ENDDO !ms
3669  ENDDO !mh
3670 
3671  ENDIF
3672 
3673  DEALLOCATE(matrix_k_i)
3674  DEALLOCATE(vm)
3675 
3678  !Store all these in equations matrices/somewhere else?????
3679 
3680  inside=.true.
3681  between=.false.
3682 
3683  source_field=>equations%INTERPOLATION%SOURCE_FIELD
3684  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
3685  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
3686  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
3687  equations_matrices=>equations%EQUATIONS_MATRICES
3688  linear_matrices=>equations_matrices%LINEAR_MATRICES
3689  equations_matrix=>linear_matrices%MATRICES(1)%PTR
3690  rhs_vector=>equations_matrices%RHS_VECTOR
3691  source_vector=>equations_matrices%SOURCE_VECTOR
3692  equations_mapping=>equations%EQUATIONS_MAPPING
3693  linear_mapping=>equations_mapping%LINEAR_MAPPING
3694  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
3695  field_var_type=field_variable%VARIABLE_TYPE
3696  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3697  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3698  source_basis=>source_field%DECOMPOSITION%DOMAIN(source_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3699  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3700  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
3701  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3702  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
3703  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3704  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3705  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3706  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3707  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3708  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
3709  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
3710  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
3711  IF(equations_set%SPECIFICATION(3)==equations_set_ale_pressure_poisson_subtype) THEN
3712  independent_field=>equations%INTERPOLATION%INDEPENDENT_FIELD
3713  independent_basis=>independent_field%DECOMPOSITION%DOMAIN(independent_field% &
3714  & decomposition%MESH_COMPONENT_NUMBER)%PTR% &
3715  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3716  CALL field_interpolation_parameters_element_get(field_mesh_velocity_set_type,element_number,equations%INTERPOLATION% &
3717  & independent_interp_parameters(field_var_type)%PTR,err,error,*999)
3718  ENDIF
3719 
3720  diff_coeff1=1.0_dp
3721  diff_coeff2=1.0_dp
3722 
3723  !Determine inside outside nodes/elements
3724  CALL field_parameter_set_data_get(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
3725  & field_input_label_set_type,input_label,err,error,*999)
3726  DO i=1,geometric_field%decomposition%domain(1)%ptr%topology%elements%maximum_number_of_element_parameters
3727  element_node_identity=geometric_field%decomposition%domain(1)%ptr%topology% &
3728  & elements%elements(element_number)%element_nodes(i)
3729  !this needs to be changed even at the right INPUT_LABEL
3730  IF(input_label(element_node_identity)<0.5_dp) THEN
3731  !labelling if !!!ELEMENT!!! is inside or outside
3732  !so even if only one node is outside node, the whole element is outside element.
3733  inside=.false.
3734  ENDIF
3735  ENDDO
3736  DO i=1,geometric_field%decomposition%domain(1)%ptr%topology%elements%maximum_number_of_element_parameters
3737  element_node_identity=geometric_field%decomposition%domain(1)%ptr%topology% &
3738  & elements%elements(element_number)%element_nodes(i)
3739  !this needs to be changed even at the right INPUT_LABEL
3740  IF(.NOT.inside) THEN
3741  IF(input_label(element_node_identity)>0.5_dp) THEN
3742  !labelling if !!!ELEMENT!!! is inside or outside
3743  !so even if only one node is outside node, the whole element is outside element.
3744  between=.true.
3745  ENDIF
3746  ENDIF
3747  ENDDO
3748 
3749  IF(inside) THEN
3750  diff_coeff1=1.0_dp
3751  diff_coeff2=1.0_dp
3752  ELSE IF(between) THEN
3753 ! set to "zero" if "between" nodes should be outside!
3754 ! DIFF_COEFF1=0.0_DP
3755  diff_coeff1=1.0_dp/1000.0_dp
3756  diff_coeff2=0.0_dp
3757  ELSE
3758  diff_coeff1=0.0_dp
3759  diff_coeff2=0.0_dp
3760  ENDIF
3761  !Loop over gauss points
3762  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
3763  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3764  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
3765  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3766  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
3767  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3768  & dependent_interp_point(field_var_type)%PTR,err,error,*999)
3769  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
3770  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
3771  IF(equations_set%SPECIFICATION(3)==equations_set_ale_pressure_poisson_subtype) THEN
3772  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3773  & independent_interp_point(field_u_variable_type)%PTR,err,error,*999)
3774  w_value(1)=equations%INTERPOLATION%INDEPENDENT_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
3775  w_value(2)=equations%INTERPOLATION%INDEPENDENT_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,no_part_deriv)
3776  IF(independent_basis%NUMBER_OF_XI==3) THEN
3777  w_value(3)=equations%INTERPOLATION%INDEPENDENT_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,no_part_deriv)
3778  END IF
3779  ELSE
3780  w_value=0.0_dp
3781  END IF
3782 
3783 ! ! ! W_VALUE=0.0_DP
3784 
3785  u_value=0.0_dp
3786  u_old=0.0_dp
3787  dxi_dx=0.0_dp
3788  dxi_dx_dx=0.0_dp
3789  delta_t=1.0_dp !need to be dependent on example file
3790  u_deriv=0.0_dp
3791  u_second=0.0_dp
3792  IF(source_vector%UPDATE_VECTOR) THEN
3793  CALL field_interpolation_parameters_element_get(field_input_data1_set_type,element_number, &
3794  & equations%INTERPOLATION%SOURCE_INTERP_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
3795  CALL field_interpolate_gauss(second_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3796  & source_interp_point(field_u_variable_type)%PTR,err,error,*999)
3797  u_value(1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
3798  u_value(2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,no_part_deriv)
3799  u_deriv(1,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s1)
3800  u_deriv(1,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s2)
3801  u_deriv(2,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s1)
3802  u_deriv(2,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s2)
3803  u_second(1,1,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s1_s1)
3804  u_second(1,1,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s1_s2)
3805  u_second(1,2,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s1_s2)
3806  u_second(1,2,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s2_s2)
3807  u_second(2,1,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s1_s1)
3808  u_second(2,1,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s1_s2)
3809  u_second(2,2,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s1_s2)
3810  u_second(2,2,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s2_s2)
3811  IF(dependent_basis%NUMBER_OF_XI==3) THEN
3812  u_value(3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,no_part_deriv)
3813  u_deriv(1,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s3)
3814  u_deriv(2,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s3)
3815  u_deriv(3,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s1)
3816  u_deriv(3,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s2)
3817  u_deriv(3,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s3)
3818  u_second(1,1,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s1_s3)
3819  u_second(1,2,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s2_s3)
3820  u_second(1,3,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s1_s3)
3821  u_second(1,3,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s2_s3)
3822  u_second(1,3,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,part_deriv_s3_s3)
3823  u_second(2,1,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s1_s3)
3824  u_second(2,2,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s2_s3)
3825  u_second(2,3,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s1_s3)
3826  u_second(2,3,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s2_s3)
3827  u_second(2,3,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,part_deriv_s3_s3)
3828  u_second(3,1,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s1_s1)
3829  u_second(3,1,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s1_s2)
3830  u_second(3,1,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s1_s3)
3831  u_second(3,2,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s1_s2)
3832  u_second(3,2,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s2_s2)
3833  u_second(3,2,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s2_s3)
3834  u_second(3,3,1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s1_s3)
3835  u_second(3,3,2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s2_s3)
3836  u_second(3,3,3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,part_deriv_s3_s3)
3837  ENDIF
3838  CALL field_interpolation_parameters_element_get(field_input_data2_set_type,element_number, &
3839  & equations%INTERPOLATION%SOURCE_INTERP_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
3840  CALL field_interpolate_gauss(second_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
3841  & source_interp_point(field_u_variable_type)%PTR,err,error,*999)
3842  u_old(1)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
3843  u_old(2)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,no_part_deriv)
3844  p_deriv(1)=equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR%VALUES(1,part_deriv_s1)
3845  p_deriv(2)=equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR%VALUES(1,part_deriv_s2)
3846  IF(dependent_basis%NUMBER_OF_XI==3) THEN
3847  u_old(3)=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,no_part_deriv)
3848  p_deriv(3)=equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR%VALUES(1,part_deriv_s3)
3849  ENDIF
3850  ENDIF
3851  !Define RHO_PARAM, density=1
3852  mu_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
3853  rho_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,no_part_deriv)
3854  !Calculate RWG.
3855 !!TODO: Think about symmetric problems.
3856  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
3857  & quadrature_scheme%GAUSS_WEIGHTS(ng)
3858  !Loop over field components
3859  mhs=0
3860  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3861  !Loop over element rows
3862 !!TODO: CHANGE ELEMENT CALCULATE TO WORK OF ns ???
3863  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
3864  mhs=mhs+1
3865  nhs=0
3866  IF(equations_matrix%UPDATE_MATRIX) THEN
3867  !Loop over element columns
3868  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
3869  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
3870  nhs=nhs+1
3871  phins=quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
3872  DO ni=1,dependent_basis%NUMBER_OF_XI
3873  pgmsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
3874  pgnsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
3875  DO mi=1,dependent_basis%NUMBER_OF_XI
3876  dxi_dx(mi,ni)=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR% &
3877  & dxi_dx(mi,ni)
3878  END DO
3879  ENDDO !ni
3880  sum=0.0_dp
3881 
3882  DO mi=1,dependent_basis%NUMBER_OF_XI
3883  DO ni=1,dependent_basis%NUMBER_OF_XI
3884  sum=sum+pgmsi(mi)*pgnsi(ni)*equations%INTERPOLATION% &
3885  & geometric_interp_point_metrics(field_u_variable_type)%PTR%GU(mi,ni)*diff_coeff1
3886  ENDDO !ni
3887  ENDDO !mi
3888  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg
3889  ENDDO !ns
3890  ENDDO !nh
3891  ENDIF
3892  IF(equations_set%SPECIFICATION(3)==equations_set_linear_pressure_poisson_subtype) THEN
3893  sum=0.0_dp
3894  sum2=0.0_dp
3895  b=0.0_dp
3896  IF(source_vector%UPDATE_VECTOR) THEN
3897  DO k=1,3
3898  !this is the temporal term (should be zero for static)
3899  IF(inside) THEN
3900  b(k)=-rho_param*((u_value(k)-u_old(k))/delta_t)
3901 ! ! this is the first viscous term
3902 ! DO K=1,3
3903 ! DO L=1,3
3904 ! B(I)=B(I)+U_DERIV(I,L)*DXI_DX_DX(L,K)
3905 ! ENDDO
3906 ! ENDDO
3907 ! this is the second viscous term
3908  DO l=1,3
3909  DO i=1,3
3910  ! this is the ALE term
3911  b(k)=b(k)+rho_param*w_value(i)*u_deriv(k,l)*dxi_dx(l,i)
3912  DO h=1,3
3913  b(k)=b(k)+mu_param*u_second(k,l,h)*dxi_dx(l,i)*dxi_dx(h,i)
3914  ENDDO
3915  ENDDO
3916  ENDDO
3917  ELSE
3918  b(k)=0.0_dp
3919  ENDIF
3920  !eventually it gets combined to the RHS (source) term
3921  DO j=1,3
3922  sum=sum+b(k)*pgmsi(j)*dxi_dx(j,k)
3923  ENDDO
3924  ENDDO
3925  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+sum*rwg
3926  ENDIF
3927  ELSE IF(equations_set%SPECIFICATION(3)==equations_set_nonlinear_pressure_poisson_subtype.OR. &
3928  & equations_set%SPECIFICATION(3)==equations_set_ale_pressure_poisson_subtype.OR. &
3929  & equations_set%SPECIFICATION(3)==equations_set_fitted_pressure_poisson_subtype) THEN
3930  sum=0.0_dp
3931  sum2=0.0_dp
3932  b=0.0_dp
3933  IF(source_vector%UPDATE_VECTOR) THEN
3934  DO k=1,3
3935  IF(inside) THEN
3936  !this is the temporal term (should be zero for static)
3937  b(k)=-rho_param*((u_value(k)-u_old(k))/delta_t)
3938  !this is the nonlinear term (should be zero for Stokes)
3939  DO l=1,3
3940  DO i=1,3
3941  !this is the nonlinear/ALE trem
3942  b(k)=b(k)-rho_param*(u_value(i)-w_value(i))*u_deriv(k,l)*dxi_dx(l,i)
3943  DO h=1,3
3944  !this is the viscous terms
3945  b(k)=b(k)+mu_param*u_second(k,l,h)*dxi_dx(l,i)*dxi_dx(h,i)
3946  ENDDO
3947  ENDDO
3948  ENDDO
3949  ELSE
3950  b(k)=0.0_dp
3951  ENDIF
3952  !now bring the test function in
3953  DO j=1,3
3954  sum=sum+b(k)*pgmsi(j)*dxi_dx(j,k)*diff_coeff2
3955  ENDDO
3956  ENDDO
3957  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
3958  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_pressure_poisson_three_dim_1) THEN
3959  x(1) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,1)
3960  x(2) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,1)
3961  x(3) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,1)
3962  sum2=-4.0_dp*pi*pi/100.0*(3.0_dp*sin(2.0_dp*pi*x(1)/10.0_dp)*sin(2.0_dp*pi*x(2)/10.0)* &
3963  & sin(2.0_dp*pi*x(3)/10.0_dp)-6.0_dp*rho_param*cos(2.0_dp*pi*x(1)/10.0_dp)**2+ &
3964  & 8.0_dp*rho_param*cos(2.0_dp*pi*x(1)/10.0_dp)**2*cos(2.0_dp*pi*x(3)/10.0_dp)**2- &
3965  & 2.0_dp*rho_param*cos(2.0_dp*pi*x(3)/10.0_dp)**2+2.0_dp*rho_param*cos(2.0_dp*pi*x(1)/10.0_dp)**2* &
3966  & cos(2.0_dp*pi*x(2)/10.0_dp)**2+4.0_dp*rho_param*cos(2.0_dp*pi*x(2)/10.0_dp)**2- &
3967  & 2.0_dp*rho_param*cos(2.0_dp*pi*x(2)/10.0_dp)**2*cos(2.0_dp*pi*x(3)/10.0_dp)**2)
3968  ELSE IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_pressure_poisson_three_dim_2) THEN
3969  x(1) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,1)
3970  x(2) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(2,1)
3971  x(3) = equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR%VALUES(3,1)
3972  sum2=-12.0_dp* sin(2.0_dp*pi*x(1)/10.0_dp)*pi*pi/100.0_dp*sin(2.0_dp*pi*x(2)/10.0_dp)* &
3973  & sin(2.0_dp*pi*x(3)/10.0_dp)
3974  sum=0.0_dp
3975  ELSE
3976  CALL flagerror("Not implemented.",err,error,*999)
3977  ENDIF
3978  ENDIF
3979  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+sum*rwg &
3980  & -sum2*quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*rwg
3981  ENDIF
3982  ELSE
3983  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
3984  ENDIF
3985  ENDDO !ms
3986  ENDDO !mh
3987  ENDDO !ng
3988  !Scale factor adjustment
3989  IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
3990  CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
3991  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
3992  mhs=0
3993  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
3994  !Loop over element rows
3995  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
3996  mhs=mhs+1
3997  nhs=0
3998  IF(equations_matrix%UPDATE_MATRIX) THEN
3999  !Loop over element columns
4000  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
4001  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4002  nhs=nhs+1
4003  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
4004  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
4005  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
4006  ENDDO !ns
4007  ENDDO !nh
4008  ENDIF
4009  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
4010  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
4011  IF(source_vector%UPDATE_VECTOR) source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector% &
4012  & element_vector%VECTOR(mhs)*equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)% &
4013  & ptr%SCALE_FACTORS(ms,mh)
4014  ENDDO !ms
4015  ENDDO !mh
4016  ENDIF
4018  CALL flagerror("Not implemented.",err,error,*999)
4020  CALL flagerror("Can not calculate finite element stiffness matrices for a nonlinear source.",err,error,*999)
4022  CALL flagerror("Can not calculate finite element stiffness matrices for a nonlinear source.",err,error,*999)
4023  CASE DEFAULT
4024  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
4025  & " is not valid for a Poisson equation type of a classical field equations set class."
4026  CALL flagerror(local_error,err,error,*999)
4027  END SELECT
4028  ELSE
4029  CALL flagerror("Equations set equations is not associated.",err,error,*999)
4030  ENDIF
4031  ELSE
4032  CALL flagerror("Equations set is not associated.",err,error,*999)
4033  ENDIF
4034 
4035  exits("POISSON_EQUATION_FINITE_ELEMENT_CALCULATE")
4036  RETURN
4037 999 errorsexits("POISSON_EQUATION_FINITE_ELEMENT_CALCULATE",err,error)
4038  RETURN 1
4040 
4041  !
4042  !================================================================================================================================
4043  !
4044 
4046  SUBROUTINE poisson_finiteelementjacobianevaluate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
4048  !Argument variables
4049  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4050  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
4051  INTEGER(INTG), INTENT(OUT) :: ERR
4052  TYPE(varying_string), INTENT(OUT) :: ERROR
4053  !Local Variables
4054  INTEGER(INTG) FIELD_VAR_TYPE,ng,mh,mhs,ms,nh,nhs,ns
4055  REAL(DP) :: B_PARAM,C_PARAM,RWG,U_VALUE,VALUE
4056  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
4057  TYPE(equations_type), POINTER :: EQUATIONS
4058  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
4059  TYPE(equations_mapping_nonlinear_type), POINTER :: NONLINEAR_MAPPING
4060  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
4061  TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
4062  TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
4063  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
4064  TYPE(field_variable_type), POINTER :: DEPENDENT_VARIABLE,GEOMETRIC_VARIABLE
4065  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
4066  TYPE(varying_string) :: LOCAL_ERROR
4067 
4068  enters("Poisson_FiniteElementJacobianEvaluate",err,error,*999)
4069 
4070  IF(ASSOCIATED(equations_set)) THEN
4071  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
4072  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
4073  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
4074  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
4075  & err,error,*999)
4076  END IF
4077  equations=>equations_set%EQUATIONS
4078  IF(ASSOCIATED(equations)) THEN
4079  SELECT CASE(equations_set%SPECIFICATION(3))
4081  CALL flagerror("Can not evaluate a Jacobian for a Poisson equation with a linear source.",err,error,*999)
4083  CALL flagerror("Can not evaluate a Jacobian for a Poisson equation with a linear source.",err,error,*999)
4085  CALL flagerror("Can not evaluate a Jacobian for a Poisson equation with a linear source.",err,error,*999)
4087  equations_matrices=>equations%EQUATIONS_MATRICES
4088  nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
4089  jacobian_matrix=>nonlinear_matrices%JACOBIANS(1)%PTR
4090  IF(jacobian_matrix%UPDATE_JACOBIAN) THEN
4091  !Store all these in equations matrices/somewhere else?????
4092  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
4093  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
4094  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
4095  equations_mapping=>equations%EQUATIONS_MAPPING
4096  nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
4097  dependent_variable=>nonlinear_mapping%RESIDUAL_VARIABLES(1)%PTR
4098  field_var_type=dependent_variable%VARIABLE_TYPE
4099  geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
4100  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4101  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4102  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4103  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4104  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
4105  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4106  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
4107  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4108  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4109  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4110  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4111  !Loop over gauss points
4112  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
4113  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4114  & dependent_interp_point(field_var_type)%PTR,err,error,*999)
4115  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4116  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
4117  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
4118  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
4119  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4120  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
4121  !Calculate RWG.
4122  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
4123  & quadrature_scheme%GAUSS_WEIGHTS(ng)
4124  !Find material parameters and u value at this Gauss point
4125  c_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4126  & values(geometric_variable%NUMBER_OF_COMPONENTS+3,no_part_deriv)
4127  u_value=equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR%VALUES(1,no_part_deriv)
4128  !Loop over field components
4129  mhs=0
4130  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4131  !Loop over element rows
4132  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4133  mhs=mhs+1
4134  nhs=0
4135  !Loop over element columns
4136  DO nh=1,dependent_variable%NUMBER_OF_COMPONENTS
4137  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4138  nhs=nhs+1
4139  VALUE=-2.0_dp*c_param*quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
4140  & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*u_value
4141  jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+VALUE*rwg
4142  ENDDO !ns
4143  ENDDO !nh
4144  ENDDO !ms
4145  ENDDO !mh
4146  ENDDO !ng
4147  ENDIF
4149  equations_matrices=>equations%EQUATIONS_MATRICES
4150  nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
4151  jacobian_matrix=>nonlinear_matrices%JACOBIANS(1)%PTR
4152  IF(jacobian_matrix%UPDATE_JACOBIAN) THEN
4153  !Store all these in equations matrices/somewhere else?????
4154  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
4155  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
4156  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
4157  equations_mapping=>equations%EQUATIONS_MAPPING
4158  nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
4159  dependent_variable=>nonlinear_mapping%RESIDUAL_VARIABLES(1)%PTR
4160  field_var_type=dependent_variable%VARIABLE_TYPE
4161  geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
4162  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4163  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4164  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4165  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4166  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
4167  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4168  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
4169  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4170  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4171  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4172  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4173  !Loop over gauss points
4174  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
4175  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4176  & dependent_interp_point(field_var_type)%PTR,err,error,*999)
4177  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4178  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
4179  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
4180  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
4181  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4182  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
4183  !Calculate RWG.
4184  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
4185  & quadrature_scheme%GAUSS_WEIGHTS(ng)
4186  !Find material parameter and u value at this Gauss point
4187  b_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4188  & values(geometric_variable%NUMBER_OF_COMPONENTS+2,no_part_deriv)
4189  c_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4190  & values(geometric_variable%NUMBER_OF_COMPONENTS+3,no_part_deriv)
4191  u_value=equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR%VALUES(1,no_part_deriv)
4192  !Loop over field components
4193  mhs=0
4194  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4195  !Loop over element rows
4196  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4197  mhs=mhs+1
4198  nhs=0
4199  !Loop over element columns
4200  DO nh=1,dependent_variable%NUMBER_OF_COMPONENTS
4201  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4202  nhs=nhs+1
4203  VALUE=-b_param*c_param*quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
4204  & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*exp(c_param*u_value)
4205  jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+VALUE*rwg
4206  ENDDO !ns
4207  ENDDO !nh
4208  ENDDO !ms
4209  ENDDO !mh
4210  ENDDO !ng
4211  ENDIF
4212  CASE DEFAULT
4213  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
4214  & " is not valid for a Poisson equation type of a classical field equations set class."
4215  CALL flagerror(local_error,err,error,*999)
4216  END SELECT
4217  ELSE
4218  CALL flagerror("Equations set equations is not associated.",err,error,*999)
4219  ENDIF
4220  ELSE
4221  CALL flagerror("Equations set is not associated.",err,error,*999)
4222  ENDIF
4223 
4224  exits("Poisson_FiniteElementJacobianEvaluate")
4225  RETURN
4226 999 errors("Poisson_FiniteElementJacobianEvaluate",err,error)
4227  exits("Poisson_FiniteElementJacobianEvaluate")
4228  RETURN 1
4229 
4231 
4232  !
4233  !================================================================================================================================
4234  !
4235 
4237  SUBROUTINE poisson_finiteelementresidualevaluate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
4239  !Argument variables
4240  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4241  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
4242  INTEGER(INTG), INTENT(OUT) :: ERR
4243  TYPE(varying_string), INTENT(OUT) :: ERROR
4244  !Local Variables
4245  INTEGER(INTG) FIELD_VAR_TYPE,ng,mh,mhs,ms,nj,nh,nhs,ni,ns
4246  REAL(DP) :: A_PARAM,B_PARAM,C_PARAM,K_PARAM,RWG,SUM1,SUM2,PGMJ(3),PGNJ(3),U_VALUE,WG
4247  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
4248  TYPE(equations_type), POINTER :: EQUATIONS
4249  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
4250  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
4251  TYPE(equations_mapping_nonlinear_type), POINTER :: NONLINEAR_MAPPING
4252  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
4253  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
4254  TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
4255  TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
4256  TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
4257  TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
4258  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
4259  TYPE(field_variable_type), POINTER :: DEPENDENT_VARIABLE,GEOMETRIC_VARIABLE
4260  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
4261  TYPE(varying_string) :: LOCAL_ERROR
4262 
4263  enters("Poisson_FiniteElementResidualEvaluate",err,error,*999)
4264 
4265  IF(ASSOCIATED(equations_set)) THEN
4266  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
4267  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
4268  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
4269  CALL flagerror("Equations set specification must have three entries for a Poisson type equations set.", &
4270  & err,error,*999)
4271  END IF
4272  equations=>equations_set%EQUATIONS
4273  IF(ASSOCIATED(equations)) THEN
4274  SELECT CASE(equations_set%SPECIFICATION(3))
4276  CALL flagerror("Can not evaluate a residual for a Poisson equation with a linear source.",err,error,*999)
4278  CALL flagerror("Can not evaluate a residual for a Poisson equation with a linear source.",err,error,*999)
4280  CALL flagerror("Can not evaluate a residual for a Poisson equation with a linear source.",err,error,*999)
4282  !Store all these in equations matrices/somewhere else?????
4283  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
4284  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
4285  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
4286  equations_matrices=>equations%EQUATIONS_MATRICES
4287  linear_matrices=>equations_matrices%LINEAR_MATRICES
4288  nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
4289  equations_matrix=>linear_matrices%MATRICES(1)%PTR
4290  rhs_vector=>equations_matrices%RHS_VECTOR
4291  source_vector=>equations_matrices%SOURCE_VECTOR
4292  equations_mapping=>equations%EQUATIONS_MAPPING
4293  linear_mapping=>equations_mapping%LINEAR_MAPPING
4294  nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
4295  dependent_variable=>nonlinear_mapping%RESIDUAL_VARIABLES(1)%PTR
4296  field_var_type=dependent_variable%VARIABLE_TYPE
4297  geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
4298  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4299  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4300  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4301  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4302  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
4303  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4304  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
4305  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4306  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4307  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4308  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4309  !Loop over gauss points
4310  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
4311  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4312  & dependent_interp_point(field_var_type)%PTR,err,error,*999)
4313  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4314  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
4315  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
4316  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
4317  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4318  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
4319  !Calculate RWG.
4320 !!TODO: Think about symmetric problems.
4321  wg=quadrature_scheme%GAUSS_WEIGHTS(ng)
4322  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN*wg
4323  !Loop over field components
4324  IF(equations_matrix%FIRST_ASSEMBLY) THEN
4325  IF(equations_matrix%UPDATE_MATRIX) THEN
4326  b_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4327  & values(geometric_variable%NUMBER_OF_COMPONENTS+2,no_part_deriv)
4328  mhs=0
4329  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4330  !Loop over element rows
4331  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4332  mhs=mhs+1
4333  nhs=0
4334  !Loop over element columns
4335  DO nh=1,dependent_variable%NUMBER_OF_COMPONENTS
4336  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4337  nhs=nhs+1
4338  sum1=0.0_dp
4339  DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
4340  pgmj(nj)=0.0_dp
4341  pgnj(nj)=0.0_dp
4342  DO ni=1,geometric_basis%NUMBER_OF_XI
4343  pgmj(nj)=pgmj(nj)+ &
4344  & quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)* &
4345  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4346  pgnj(nj)=pgnj(nj)+ &
4347  & quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)* &
4348  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4349  ENDDO !ni
4350  k_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(nj,no_part_deriv)
4351  sum1=sum1+k_param*pgmj(nj)*pgnj(nj)
4352  ENDDO !nj
4353  sum2=b_param*quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
4354  & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
4355  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+ &
4356  & (sum1+sum2)*rwg
4357  ENDDO !ns
4358  ENDDO !nh
4359  ENDDO !ms
4360  ENDDO !mh
4361  ENDIF
4362  ENDIF
4363  IF(rhs_vector%FIRST_ASSEMBLY) THEN
4364  IF(rhs_vector%UPDATE_VECTOR) THEN
4365  a_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4366  & values(geometric_variable%NUMBER_OF_COMPONENTS+1,no_part_deriv)
4367  mhs=0
4368  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4369  !Loop over element rows
4370  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4371  mhs=mhs+1
4372  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)+ &
4373  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*a_param*rwg
4374  ENDDO !ms
4375  ENDDO !mh
4376  ENDIF
4377  ENDIF
4378  IF(nonlinear_matrices%UPDATE_RESIDUAL) THEN
4379  c_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4380  & values(geometric_variable%NUMBER_OF_COMPONENTS+3,no_part_deriv)
4381  u_value=equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_u_variable_type)%PTR% &
4382  & values(1,no_part_deriv)
4383  mhs=0
4384  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4385  !Loop over element rows
4386  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4387  mhs=mhs+1
4388  nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)- &
4389  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*c_param*u_value**2*rwg
4390  ENDDO !ms
4391  ENDDO !mh
4392  ENDIF
4393  ENDDO !ng
4395  !Store all these in equations matrices/somewhere else?????
4396  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
4397  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
4398  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
4399  equations_matrices=>equations%EQUATIONS_MATRICES
4400  linear_matrices=>equations_matrices%LINEAR_MATRICES
4401  nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
4402  equations_matrix=>linear_matrices%MATRICES(1)%PTR
4403  rhs_vector=>equations_matrices%RHS_VECTOR
4404  source_vector=>equations_matrices%SOURCE_VECTOR
4405  equations_mapping=>equations%EQUATIONS_MAPPING
4406  linear_mapping=>equations_mapping%LINEAR_MAPPING
4407  nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
4408  dependent_variable=>nonlinear_mapping%RESIDUAL_VARIABLES(1)%PTR
4409  field_var_type=dependent_variable%VARIABLE_TYPE
4410  geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
4411  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4412  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4413  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4414  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4415  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
4416  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4417  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
4418  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4419  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4420  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4421  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4422  !Loop over gauss points
4423  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
4424  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4425  & dependent_interp_point(field_var_type)%PTR,err,error,*999)
4426  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4427  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
4428  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
4429  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
4430  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4431  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
4432  !Calculate RWG.
4433 !!TODO: Think about symmetric problems.
4434  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
4435  & quadrature_scheme%GAUSS_WEIGHTS(ng)
4436  !Loop over field components
4437  IF(equations_matrix%FIRST_ASSEMBLY) THEN
4438  IF(equations_matrix%UPDATE_MATRIX) THEN
4439  mhs=0
4440  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4441  !Loop over element rows
4442  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4443  mhs=mhs+1
4444  nhs=0
4445  !Loop over element columns
4446  DO nh=1,dependent_variable%NUMBER_OF_COMPONENTS
4447  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4448  nhs=nhs+1
4449  sum1=0.0_dp
4450  DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
4451  pgmj(nj)=0.0_dp
4452  pgnj(nj)=0.0_dp
4453  DO ni=1,geometric_basis%NUMBER_OF_XI
4454  pgmj(nj)=pgmj(nj)+ &
4455  & quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)* &
4456  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4457  pgnj(nj)=pgnj(nj)+ &
4458  & quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)* &
4459  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4460  ENDDO !ni
4461  k_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(nj,no_part_deriv)
4462  sum1=sum1+k_param*pgmj(nj)*pgnj(nj)
4463  ENDDO !nj
4464  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum1*rwg
4465  ENDDO !ns
4466  ENDDO !nh
4467  ENDDO !ms
4468  ENDDO !mh
4469  ENDIF
4470  ENDIF
4471  IF(rhs_vector%FIRST_ASSEMBLY) THEN
4472  IF(rhs_vector%UPDATE_VECTOR) THEN
4473  a_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4474  & values(geometric_variable%NUMBER_OF_COMPONENTS+1,no_part_deriv)
4475  mhs=0
4476  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4477  !Loop over element rows
4478  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4479  mhs=mhs+1
4480  rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)+ &
4481  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*a_param*rwg
4482  ENDDO !ms
4483  ENDDO !mh
4484  ENDIF
4485  ENDIF
4486  IF(nonlinear_matrices%UPDATE_RESIDUAL) THEN
4487  b_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4488  & values(geometric_variable%NUMBER_OF_COMPONENTS+2,no_part_deriv)
4489  c_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4490  & values(geometric_variable%NUMBER_OF_COMPONENTS+3,no_part_deriv)
4491  u_value=equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_u_variable_type)%PTR% &
4492  & values(1,no_part_deriv)
4493 !!TODO: Handle floating point exceptions better
4494  IF((c_param*u_value)>20000.0_dp) THEN
4495  local_error="The value of "//trim(number_to_vstring(c_param*u_value,"*",err,error))// &
4496  & " is out of range for an exponential function."
4497  CALL flagerror(local_error,err,error,*999)
4498  ENDIF
4499  mhs=0
4500  DO mh=1,dependent_variable%NUMBER_OF_COMPONENTS
4501  !Loop over element rows
4502  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4503  mhs=mhs+1
4504  nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)- &
4505  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*b_param*exp(c_param*u_value)*rwg
4506  ENDDO !ms
4507  ENDDO !mh
4508  ENDIF
4509  ENDDO !ng
4510  CASE DEFAULT
4511  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
4512  & " is not valid for a Poisson equation type of a classical field equations set class."
4513  CALL flagerror(local_error,err,error,*999)
4514  END SELECT
4515  ELSE
4516  CALL flagerror("Equations set equations is not associated.",err,error,*999)
4517  ENDIF
4518  ELSE
4519  CALL flagerror("Equations set is not associated.",err,error,*999)
4520  ENDIF
4521 
4522  exits("Poisson_FiniteElementResidualEvaluate")
4523  RETURN
4524 999 errors("Poisson_FiniteElementResidualEvaluate",err,error)
4525  exits("Poisson_FiniteElementResidualEvaluate")
4526  RETURN 1
4527 
4529 
4530  !
4531  !================================================================================================================================
4532  !
4533 
4535  SUBROUTINE poisson_problemspecificationset(problem,problemSpecification,err,error,*)
4537  !Argument variables
4538  TYPE(problem_type), POINTER :: problem
4539  INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
4540  INTEGER(INTG), INTENT(OUT) :: err
4541  TYPE(varying_string), INTENT(OUT) :: error
4542  !Local Variables
4543  TYPE(varying_string) :: localError
4544  INTEGER(INTG) :: problemSubtype
4545 
4546  enters("Poisson_ProblemSpecificationSet",err,error,*999)
4547 
4548  IF(ASSOCIATED(problem)) THEN
4549  IF(SIZE(problemspecification,1)==3) THEN
4550  problemsubtype=problemspecification(3)
4551  SELECT CASE(problemsubtype)
4559  !ok
4560  CASE DEFAULT
4561  localerror="The third problem specification of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
4562  & " is not valid for a Poisson problem."
4563  CALL flagerror(localerror,err,error,*999)
4564  END SELECT
4565  IF(ALLOCATED(problem%specification)) THEN
4566  CALL flagerror("Problem specification is already allocated.",err,error,*999)
4567  ELSE
4568  ALLOCATE(problem%specification(3),stat=err)
4569  IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
4570  END IF
4571  problem%specification(1:3)=[problem_classical_field_class,problem_poisson_equation_type,problemsubtype]
4572  ELSE
4573  CALL flagerror("Poisson problem specification must have three entries.",err,error,*999)
4574  END IF
4575  ELSE
4576  CALL flagerror("Problem is not associated.",err,error,*999)
4577  END IF
4578 
4579  exits("Poisson_ProblemSpecificationSet")
4580  RETURN
4581 999 errorsexits("Poisson_ProblemSpecificationSet",err,error)
4582  RETURN 1
4583 
4584  END SUBROUTINE poisson_problemspecificationset
4585 
4586  !
4587  !================================================================================================================================
4588  !
4589 
4591  SUBROUTINE poisson_problemextracellularbidomainsetup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
4593  !Argument variables
4594  TYPE(problem_type), POINTER :: PROBLEM
4595  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
4596  INTEGER(INTG), INTENT(OUT) :: ERR
4597  TYPE(varying_string), INTENT(OUT) :: ERROR
4598  !Local Variables
4599  TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
4600  TYPE(solver_type), POINTER :: SOLVER
4601  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4602  TYPE(solvers_type), POINTER :: SOLVERS
4603  TYPE(varying_string) :: LOCAL_ERROR
4604 
4605  enters("Poisson_ProblemExtracellularBidomainSetup",err,error,*999)
4606 
4607  NULLIFY(control_loop)
4608  NULLIFY(solver)
4609  NULLIFY(solver_equations)
4610  NULLIFY(solvers)
4611  IF(ASSOCIATED(problem)) THEN
4612  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
4613  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4614  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
4615  CALL flagerror("Problem specification must have three entries for a Poisson problem.",err,error,*999)
4616  END IF
4617  IF(problem%SPECIFICATION(3)==problem_extracellular_bidomain_poisson_subtype) THEN
4618  SELECT CASE(problem_setup%SETUP_TYPE)
4620  SELECT CASE(problem_setup%ACTION_TYPE)
4622  !Do nothing????
4624  !Do nothing????
4625  CASE DEFAULT
4626  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4627  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4628  & " is invalid for a extracellular bidomain Poisson equation."
4629  CALL flagerror(local_error,err,error,*999)
4630  END SELECT
4632  SELECT CASE(problem_setup%ACTION_TYPE)
4634  !Set up a simple control loop
4635  CALL control_loop_create_start(problem,control_loop,err,error,*999)
4637  !Finish the control loops
4638  control_loop_root=>problem%CONTROL_LOOP
4639  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4640  CALL control_loop_create_finish(control_loop,err,error,*999)
4641  CASE DEFAULT
4642  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4643  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4644  & " is invalid for a extracellular bidomain Poisson equation."
4645  CALL flagerror(local_error,err,error,*999)
4646  END SELECT
4648  !Get the control loop
4649  control_loop_root=>problem%CONTROL_LOOP
4650  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4651  SELECT CASE(problem_setup%ACTION_TYPE)
4653  !Start the solvers creation
4654  CALL solvers_create_start(control_loop,solvers,err,error,*999)
4655  CALL solvers_number_set(solvers,1,err,error,*999)
4656  !Set the solver to be a linear solver
4657  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
4658  !Start the linear solver creation
4659  CALL solver_type_set(solver,solver_linear_type,err,error,*999)
4660  !Set solver defaults
4661  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
4663  !Get the solvers
4664  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
4665  !Finish the solvers creation
4666  CALL solvers_create_finish(solvers,err,error,*999)
4667  CASE DEFAULT
4668  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4669  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4670  & " is invalid for a extracellular bidomain Poisson equation."
4671  CALL flagerror(local_error,err,error,*999)
4672  END SELECT
4674  SELECT CASE(problem_setup%ACTION_TYPE)
4676  !Get the control loop
4677  control_loop_root=>problem%CONTROL_LOOP
4678  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4679  !Get the solver
4680  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
4681  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
4682  !Create the solver equations
4683  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
4684  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
4685  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_static,err,error,*999)
4686  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
4688  !Get the control loop
4689  control_loop_root=>problem%CONTROL_LOOP
4690  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4691  !Get the solver equations
4692  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
4693  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
4694  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
4695  !Finish the solver equations creation
4696  CALL solver_equations_create_finish(solver_equations,err,error,*999)
4697  CASE DEFAULT
4698  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4699  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4700  & " is invalid for a extracellular bidomain Poisson equation."
4701  CALL flagerror(local_error,err,error,*999)
4702  END SELECT
4703  CASE DEFAULT
4704  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4705  & " is invalid for a extracellular bidomain Poisson equation."
4706  CALL flagerror(local_error,err,error,*999)
4707  END SELECT
4708  ELSE
4709  local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
4710  & " does not equal a extracellular bidomain Poisson equation subtype."
4711  CALL flagerror(local_error,err,error,*999)
4712  ENDIF
4713  ELSE
4714  CALL flagerror("Problem is not associated.",err,error,*999)
4715  ENDIF
4716 
4717  exits("Poisson_ProblemExtracellularBidomainSetup")
4718  RETURN
4719 999 errors("Poisson_ProblemExtracellularBidomainSetup",err,error)
4720  exits("Poisson_ProblemExtracellularBidomainSetup")
4721  RETURN 1
4722 
4724 
4725  !
4726  !================================================================================================================================
4727  !
4728 
4730  SUBROUTINE poisson_problemlinearsourcesetup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
4732  !Argument variables
4733  TYPE(problem_type), POINTER :: PROBLEM
4734  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
4735  INTEGER(INTG), INTENT(OUT) :: ERR
4736  TYPE(varying_string), INTENT(OUT) :: ERROR
4737  !Local Variables
4738  TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
4739  TYPE(solver_type), POINTER :: SOLVER
4740  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4741  TYPE(solvers_type), POINTER :: SOLVERS
4742  TYPE(varying_string) :: LOCAL_ERROR
4743 
4744  enters("Poisson_ProblemLinearSourceSetup",err,error,*999)
4745 
4746  NULLIFY(control_loop)
4747  NULLIFY(solver)
4748  NULLIFY(solver_equations)
4749  NULLIFY(solvers)
4750  IF(ASSOCIATED(problem)) THEN
4751  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
4752  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4753  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
4754  CALL flagerror("Problem specification must have three entries for a Poisson problem.",err,error,*999)
4755  END IF
4756  IF(problem%SPECIFICATION(3)==problem_linear_source_poisson_subtype) THEN
4757  SELECT CASE(problem_setup%SETUP_TYPE)
4759  SELECT CASE(problem_setup%ACTION_TYPE)
4761  !Do nothing????
4763  !Do nothing????
4764  CASE DEFAULT
4765  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4766  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4767  & " is invalid for a linear source Poisson equation."
4768  CALL flagerror(local_error,err,error,*999)
4769  END SELECT
4771  SELECT CASE(problem_setup%ACTION_TYPE)
4773  !Set up a simple control loop
4774  CALL control_loop_create_start(problem,control_loop,err,error,*999)
4776  !Finish the control loops
4777  control_loop_root=>problem%CONTROL_LOOP
4778  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4779  CALL control_loop_create_finish(control_loop,err,error,*999)
4780  CASE DEFAULT
4781  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4782  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4783  & " is invalid for a linear source Poisson equation."
4784  CALL flagerror(local_error,err,error,*999)
4785  END SELECT
4787  !Get the control loop
4788  control_loop_root=>problem%CONTROL_LOOP
4789  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4790  SELECT CASE(problem_setup%ACTION_TYPE)
4792  !Start the solvers creation
4793  CALL solvers_create_start(control_loop,solvers,err,error,*999)
4794  CALL solvers_number_set(solvers,1,err,error,*999)
4795  !Set the solver to be a linear solver
4796  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
4797  !Start the linear solver creation
4798  CALL solver_type_set(solver,solver_linear_type,err,error,*999)
4799  !Set solver defaults
4800  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
4802  !Get the solvers
4803  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
4804  !Finish the solvers creation
4805  CALL solvers_create_finish(solvers,err,error,*999)
4806  CASE DEFAULT
4807  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4808  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4809  & " is invalid for a linear source Poisson equation."
4810  CALL flagerror(local_error,err,error,*999)
4811  END SELECT
4813  SELECT CASE(problem_setup%ACTION_TYPE)
4815  !Get the control loop
4816  control_loop_root=>problem%CONTROL_LOOP
4817  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4818  !Get the solver
4819  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
4820  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
4821  !Create the solver equations
4822  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
4823  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
4824  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_static,err,error,*999)
4825  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
4827  !Get the control loop
4828  control_loop_root=>problem%CONTROL_LOOP
4829  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4830  !Get the solver equations
4831  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
4832  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
4833  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
4834  !Finish the solver equations creation
4835  CALL solver_equations_create_finish(solver_equations,err,error,*999)
4836  CASE DEFAULT
4837  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4838  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4839  & " is invalid for a linear source Poisson equation."
4840  CALL flagerror(local_error,err,error,*999)
4841  END SELECT
4842  CASE DEFAULT
4843  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4844  & " is invalid for a linear source Poisson equation."
4845  CALL flagerror(local_error,err,error,*999)
4846  END SELECT
4847  ELSE
4848  local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
4849  & " does not equal a linear source Poisson equation subtype."
4850  CALL flagerror(local_error,err,error,*999)
4851  ENDIF
4852  ELSE
4853  CALL flagerror("Problem is not associated.",err,error,*999)
4854  ENDIF
4855 
4856  exits("Poisson_ProblemLinearSourceSetup")
4857  RETURN
4858 999 errorsexits("Poisson_ProblemLinearSourceSetup",err,error)
4859  RETURN 1
4860 
4861  END SUBROUTINE poisson_problemlinearsourcesetup
4862 
4863 
4864  !
4865  !================================================================================================================================
4866  !
4867 
4869  SUBROUTINE poisson_problempressurepoissonsetup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
4871  !Argument variables
4872  TYPE(problem_type), POINTER :: PROBLEM
4873  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
4874  INTEGER(INTG), INTENT(OUT) :: ERR
4875  TYPE(varying_string), INTENT(OUT) :: ERROR
4876  !Local Variables
4877  TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT,SUB_LOOP
4878  TYPE(solver_type), POINTER :: SOLVER,FITTING_SOLVER
4879  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS,FITTING_SOLVER_EQUATIONS
4880  TYPE(solvers_type), POINTER :: SOLVERS
4881  TYPE(varying_string) :: LOCAL_ERROR
4882 
4883  enters("Poisson_ProblemPressurePoissonSetup",err,error,*999)
4884 
4885  NULLIFY(control_loop)
4886  NULLIFY(solver)
4887  NULLIFY(fitting_solver)
4888  NULLIFY(solver_equations)
4889  NULLIFY(fitting_solver_equations)
4890  NULLIFY(solvers)
4891  IF(ASSOCIATED(problem)) THEN
4892  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
4893  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4894  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
4895  CALL flagerror("Problem specification must have three entries for a Poisson problem.",err,error,*999)
4896  END IF
4897  IF(problem%SPECIFICATION(3)==problem_nonlinear_pressure_poisson_subtype.OR. &
4898  & problem%SPECIFICATION(3)==problem_linear_pressure_poisson_subtype.OR. &
4899  & problem%SPECIFICATION(3)==problem_ale_pressure_poisson_subtype) THEN
4900  SELECT CASE(problem_setup%SETUP_TYPE)
4902  SELECT CASE(problem_setup%ACTION_TYPE)
4904  !Do nothing????
4906  !Do nothing????
4907  CASE DEFAULT
4908  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4909  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4910  & " is invalid for a linear source Poisson equation."
4911  CALL flagerror(local_error,err,error,*999)
4912  END SELECT
4914  SELECT CASE(problem_setup%ACTION_TYPE)
4916  NULLIFY(control_loop_root)
4917  NULLIFY(sub_loop)
4918  !Set up a simple control loop with 2 subloops
4919  CALL control_loop_create_start(problem,control_loop,err,error,*999)
4920  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
4921  CALL control_loop_number_of_sub_loops_set(control_loop,1,err,error,*999)
4922  CALL control_loop_sub_loop_get(control_loop,1,sub_loop,err,error,*999)
4923  CALL control_loop_type_set(sub_loop,problem_control_while_loop_type,err,error,*999)
4924  !Set the number of iterations for the while loop to 3 for now
4925 ! WRITE(*,*)'BE CAREFUL HERE ! MAKE IT MORE GENERAL'
4926  CALL control_loop_maximum_iterations_set(sub_loop,1,err,error,*999)
4928  !Finish the control loops
4929  control_loop_root=>problem%CONTROL_LOOP
4930  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
4931  CALL control_loop_create_finish(control_loop,err,error,*999)
4932  CASE DEFAULT
4933  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
4934  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
4935  & " is invalid for a linear source Poisson equation."
4936  CALL flagerror(local_error,err,error,*999)
4937  END SELECT
4939  !Get the control loop
4940  control_loop_root=>problem%CONTROL_LOOP
4941  NULLIFY(sub_loop)
4942  NULLIFY(control_loop)
4943  CALL control_loop_sub_loop_get(control_loop_root,1,sub_loop,err,error,*999)
4944  CALL control_loop_get(sub_loop,control_loop_node,control_loop,err,error,*999)
4945 
4946  SELECT CASE(problem_setup%ACTION_TYPE)
4947  CASE(