OpenCMISS-Iron Internal API Documentation
diffusion_equation_routines.f90
Go to the documentation of this file.
1 
43 
46 
48  USE base_routines
49  USE basis_routines
51  USE constants
54  USE domain_mappings
60  USE field_routines
61  USE input_output
63  USE kinds
64  USE matrix_vector
66  USE strings
67  USE solver_routines
68  USE timer
69  USE types
70 
72 
73 #include "macros.h"
74 
75  IMPLICIT NONE
76 
77  PRIVATE
78 
79  !Module parameters
80 
81  !Module types
82 
83  !Module variables
84 
85  !Interfaces
86 
88 
90 
92 
94 
96 
98 
100 
102 
104 
106 
108 
109 !!TODO: should the following two routines really be public???
110 
112 
114 
115 
116 CONTAINS
117 
118  !
119  !================================================================================================================================
120  !
121 
122 
124  !Calculates a two-dimensional unsteady heat equation solution (diffusion coefficient is 1)
125  SUBROUTINE diffusion_boundaryconditionanalyticcalculate(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*)
127  !Argument variables
128  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
129  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
130  INTEGER(INTG), INTENT(OUT) :: ERR
131  TYPE(varying_string), INTENT(OUT) :: ERROR
132  !Local Variables
133  INTEGER(INTG) :: component_idx,deriv_idx,dim_idx,local_ny,node_idx,NUMBER_OF_DIMENSIONS,variable_idx,variable_type,version_idx
134  REAL(DP) :: VALUE,X(3),INITIAL_VALUE
135  REAL(DP), POINTER :: ANALYTIC_PARAMETERS(:),GEOMETRIC_PARAMETERS(:),MATERIALS_PARAMETERS(:)
136  TYPE(domain_type), POINTER :: DOMAIN
137  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
138  TYPE(field_type), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,EQUATIONS_SET_FIELD_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
139  TYPE(field_variable_type), POINTER :: ANALYTIC_VARIABLE,FIELD_VARIABLE,GEOMETRIC_VARIABLE,MATERIALS_VARIABLE
140  INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
141  !TYPE(VARYING_STRING) :: LOCAL_ERROR
142  INTEGER(INTG) :: GLOBAL_DERIV_INDEX,ANALYTIC_FUNCTION_TYPE,imy_matrix
143  !THESE ARE TEMPORARY VARIABLES - they need to be replace by constant field values and the current simulation time
144  REAL(DP) :: TIME,NORMAL(3),TANGENTS(3,3)
145  !CURRENT_TIME = 1.2_DP
146 
147  enters("Diffusion_BoundaryConditionAnalyticCalculate",err,error,*999)
148 
149  IF(ASSOCIATED(equations_set)) THEN
150  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
151  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
152  IF(ASSOCIATED(dependent_field)) THEN
153  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
154  IF(ASSOCIATED(geometric_field)) THEN
155  analytic_function_type=equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE
156  analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
157  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
158  NULLIFY(geometric_variable)
159  NULLIFY(geometric_parameters)
160  CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
161  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,geometric_parameters, &
162  & err,error,*999)
163  NULLIFY(analytic_variable)
164  NULLIFY(analytic_parameters)
165  IF(ASSOCIATED(analytic_field)) THEN
166  CALL field_variable_get(analytic_field,field_u_variable_type,analytic_variable,err,error,*999)
167  CALL field_parameter_set_data_get(analytic_field,field_u_variable_type,field_values_set_type, &
168  & analytic_parameters,err,error,*999)
169  ENDIF
170  NULLIFY(materials_field)
171  NULLIFY(materials_variable)
172  NULLIFY(materials_parameters)
173  IF(ASSOCIATED(equations_set%MATERIALS)) THEN
174  materials_field=>equations_set%MATERIALS%MATERIALS_FIELD
175  CALL field_variable_get(materials_field,field_u_variable_type,materials_variable,err,error,*999)
176  CALL field_parameter_set_data_get(materials_field,field_u_variable_type,field_values_set_type, &
177  & materials_parameters,err,error,*999)
178  ENDIF
179  IF(ASSOCIATED(boundary_conditions)) THEN
180  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
181  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
182  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
183  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
184  & err,error,*999)
185  END IF
186  IF(equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype)THEN
187  !If a multi-comp model, we will use the equations set field information to assign only the appropriate field
188  !variable boundary conditions
189  time=equations_set%ANALYTIC%ANALYTIC_USER_PARAMS(1)
190  !Use predetermined mapping from equations set field compartment number to field variable type
191  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
192  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
193  & field_values_set_type,equations_set_field_data,err,error,*999)
194  imy_matrix = equations_set_field_data(1)
195  DO variable_idx=0,1
196  variable_type=field_u_variable_type+(field_number_of_variable_subtypes*(imy_matrix-1))+variable_idx
197  !variable_type=DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
198  field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
199  IF(ASSOCIATED(field_variable)) THEN
200  CALL field_parameter_set_create(dependent_field,variable_type,field_analytic_values_set_type,err,error,*999)
201  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
202  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
203  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
204  IF(ASSOCIATED(domain)) THEN
205  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
206  domain_nodes=>domain%TOPOLOGY%NODES
207  IF(ASSOCIATED(domain_nodes)) THEN
208  !Loop over the local nodes excluding the ghosts.
209  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
210 !!TODO \todo We should interpolate the geometric field here and the node position.
211  DO dim_idx=1,number_of_dimensions
212  !Default to version 1 of each node derivative
213  local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
214  & nodes(node_idx)%DERIVATIVES(1)%VERSIONS(1)
215  x(dim_idx)=geometric_parameters(local_ny)
216  ENDDO !dim_idx
217  !Loop over the derivatives
218  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
219  global_deriv_index=domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX
220  CALL diffusion_analyticfunctionsevaluate(equations_set,analytic_function_type,&
221  & x,tangents,normal,time,variable_type,global_deriv_index,component_idx, &
222  & analytic_parameters,materials_parameters,VALUE,err,error,*999)
223  !Default to version 1 of each node derivative
224  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
225  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
226  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
227  & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
228  IF(mod(variable_type,field_number_of_variable_subtypes)==field_u_variable_type) THEN
229  IF(domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN
230  !If we are a boundary node then set the analytic value on the boundary
231  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field,variable_type, &
232  & local_ny,boundary_condition_fixed,VALUE,err,error,*999)
233  ELSE
234  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
235  & field_values_set_type,local_ny,VALUE,err,error,*999)
236  ENDIF
237  ENDIF
238  ENDDO !deriv_idx
239  ENDDO !node_idx
240  ELSE
241  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
242  ENDIF
243  ELSE
244  CALL flagerror("Domain topology is not associated.",err,error,*999)
245  ENDIF
246  ELSE
247  CALL flagerror("Domain is not associated.",err,error,*999)
248  ENDIF
249  ELSE
250  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
251  ENDIF
252  ENDDO !component_idx
253  CALL field_parameter_set_update_start(dependent_field,variable_type,field_analytic_values_set_type, &
254  & err,error,*999)
255  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_analytic_values_set_type, &
256  & err,error,*999)
257  CALL field_parameter_set_update_start(dependent_field,variable_type,field_values_set_type, &
258  & err,error,*999)
259  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_values_set_type, &
260  & err,error,*999)
261  ELSE
262  CALL flagerror("Field variable is not associated.",err,error,*999)
263  ENDIF
264  ENDDO !variable_idx
265  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
266  & geometric_parameters,err,error,*999)
267  ELSE
268  !for single physics diffusion problems use standard analytic calculate
269  analytic_function_type=equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE
270  time=equations_set%ANALYTIC%ANALYTIC_TIME
271  DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES
272  variable_type=dependent_field%VARIABLES(variable_idx)%VARIABLE_TYPE
273  field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
274  IF(ASSOCIATED(field_variable)) THEN
275  CALL field_parametersetensurecreated(dependent_field,variable_type,field_analytic_values_set_type, &
276  & err,error,*999)
277  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
278  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
279  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
280  IF(ASSOCIATED(domain)) THEN
281  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
282  domain_nodes=>domain%TOPOLOGY%NODES
283  IF(ASSOCIATED(domain_nodes)) THEN
284  !Loop over the local nodes excluding the ghosts.
285  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
286 !!TODO \todo We should interpolate the geometric field here and the node position.
287  DO dim_idx=1,number_of_dimensions
288  !Default to version 1 of each node derivative
289  local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
290  & nodes(node_idx)%DERIVATIVES(1)%VERSIONS(1)
291  x(dim_idx)=geometric_parameters(local_ny)
292  ENDDO !dim_idx
293  !Loop over the derivatives
294  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
295  global_deriv_index=domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX
296  CALL diffusion_analyticfunctionsevaluate(equations_set,analytic_function_type,&
297  & x,tangents,normal,0.0_dp,variable_type,global_deriv_index,component_idx, &
298  & analytic_parameters,materials_parameters,initial_value,err,error,*999)
299  CALL diffusion_analyticfunctionsevaluate(equations_set,analytic_function_type,&
300  & x,tangents,normal,time,variable_type,global_deriv_index,component_idx, &
301  & analytic_parameters,materials_parameters,VALUE,err,error,*999)
302  DO version_idx=1,domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%numberOfVersions
303  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
304  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(version_idx)
305  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
306  & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
307  IF(variable_type==field_u_variable_type) THEN
308  IF(domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN
309  !If we are a boundary node then set the analytic value on the boundary
310  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field,variable_type, &
311  & local_ny,boundary_condition_fixed,initial_value,err,error,*999)
312  ELSE
313  !Set the initial condition.
314  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
315  & field_values_set_type,local_ny,initial_value,err,error,*999)
316  ENDIF
317  ENDIF
318  ENDDO !version_idx
319  ENDDO !deriv_idx
320  ENDDO !node_idx
321  ELSE
322  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
323  ENDIF
324  ELSE
325  CALL flagerror("Domain topology is not associated.",err,error,*999)
326  ENDIF
327  ELSE
328  CALL flagerror("Domain is not associated.",err,error,*999)
329  ENDIF
330  ELSE
331  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
332  ENDIF
333  ENDDO !component_idx
334  CALL field_parameter_set_update_start(dependent_field,variable_type,field_analytic_values_set_type, &
335  & err,error,*999)
336  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_analytic_values_set_type, &
337  & err,error,*999)
338  ELSE
339  CALL flagerror("Field variable is not associated.",err,error,*999)
340  ENDIF
341  ENDDO !variable_idx
342  ENDIF
343  ENDIF
344  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
345  & geometric_parameters,err,error,*999)
346  ELSE
347  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
348  ENDIF
349  ELSE
350  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
351  ENDIF
352  ELSE
353  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
354  ENDIF
355  ELSE
356  CALL flagerror("Equations set is not associated.",err,error,*999)
357  ENDIF
358 
359  exits("Diffusion_BoundaryConditionAnalyticCalculate")
360  RETURN
361 999 errors("Diffusion_BoundaryConditionAnalyticCalculate",err,error)
362  exits("Diffusion_BoundaryConditionAnalyticCalculate")
363  RETURN 1
364 
366 
367 
368  !
369  !================================================================================================================================
370  !
372  SUBROUTINE diffusion_analyticfunctionsevaluate(EQUATIONS_SET,ANALYTIC_FUNCTION_TYPE,X, &
373  & tangents,normal,time,variable_type,global_derivative,component_number,analytic_parameters,materials_parameters, &
374  & VALUE,err,error,*)
376  !Argument variables
377  TYPE(equations_set_type), POINTER, INTENT(IN) :: EQUATIONS_SET
378  INTEGER(INTG), INTENT(IN) :: ANALYTIC_FUNCTION_TYPE
379  REAL(DP), INTENT(IN) :: X(:)
380  REAL(DP), INTENT(IN) :: TANGENTS(:,:)
381  REAL(DP), INTENT(IN) :: NORMAL(:)
382  REAL(DP), INTENT(IN) :: TIME
383  INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
384  INTEGER(INTG), INTENT(IN) :: GLOBAL_DERIVATIVE
385  INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
386  REAL(DP), INTENT(IN) :: ANALYTIC_PARAMETERS(:)
387  REAL(DP), INTENT(IN) :: MATERIALS_PARAMETERS(:)
388  REAL(DP), INTENT(OUT) :: VALUE
389  INTEGER(INTG), INTENT(OUT) :: ERR
390  TYPE(varying_string), INTENT(OUT) :: ERROR
391  !Local variables
392  REAL(DP) :: k,phi,A,B,C,D,A1,A2,A3,A4
393  REAL(DP) :: A_PARAM,B_PARAM,C_PARAM,K_PARAM,L_PARAM,CONST_PARAM,BETA_PARAM,LAMBDA_PARAM,MU_PARAM
394  INTEGER(INTG) :: EQUATIONS_SUBTYPE
395  TYPE(varying_string) :: LOCAL_ERROR
396 
397  !These are parameters for the analytical solution
398 ! k = 1.0_DP !this is a time decay constant for the exponential term
399 ! phi = 0.785398163397_DP !pi/4 - this sets the orientation of the solution relative to the axes
400  k = 10.0_dp !this is a time decay constant for the exponential term
401  phi = 1.0_dp !pi/4 - this sets the orientation of the solution relative to the axes
402 
403  !Solution parameters for
404  a1 = 0.4_dp
405  a2 = 0.3_dp
406  a3 = 0.2_dp
407  a4 = 0.1_dp
408 
409  enters("Diffusion_AnalyticFunctionsEvaluate",err,error,*999)
410 
411  IF(.NOT.ASSOCIATED(equations_set)) THEN
412  CALL flagerror("Equations set is not associated.",err,error,*999)
413  ELSE
414  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
415  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
416  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
417  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
418  & err,error,*999)
419  ELSE
420  equations_subtype=equations_set%SPECIFICATION(3)
421  END IF
422  END IF
423  SELECT CASE(equations_subtype)
425  SELECT CASE(analytic_function_type)
427  !For \del u/\del t = a \del^2 u/\del x^2
428  !u(x,t)=A.exp(-4.\mu^2.t)cos(\mu.x+B)+C
429  !see http://eqworld.ipmnet.ru/en/solutions/lpde/lpde101.pdf
430  !OpenCMISS has \del u/\del t + k \del^2 u/\del x^2 = 0, thereform with \mu=2.\pi/L we have
431  !u(x,t)=A.exp(4.\pi^2.k.t/L^2)cos(2.\pi.x/L+B)+C
432  k_param=materials_parameters(1)
433  a_param=analytic_parameters(1)
434  b_param=analytic_parameters(2)
435  c_param=analytic_parameters(3)
436  l_param=analytic_parameters(4)
437  SELECT CASE(variable_type)
438  CASE(field_u_variable_type)
439  SELECT CASE(global_derivative)
440  CASE(no_global_deriv)
441  VALUE=a_param*exp(4.0_dp*pi**2*k_param*time/l_param**2)*cos(2.0_dp*pi*x(1)/l_param+b_param)+c_param
442  CASE(global_deriv_s1)
443  CALL flagerror("Not implemented.",err,error,*999)
444  CASE DEFAULT
445  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
446  & " is invalid."
447  CALL flagerror(local_error,err,error,*999)
448  END SELECT
449  CASE(field_deludeln_variable_type)
450  SELECT CASE(global_derivative)
451  CASE(no_global_deriv)
452  VALUE=0.0_dp
453  CASE(global_deriv_s1)
454  CALL flagerror("Not implemented.",err,error,*999)
455  CASE DEFAULT
456  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
457  & " is invalid."
458  CALL flagerror(local_error,err,error,*999)
459  END SELECT
460  CASE DEFAULT
461  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
462  & " is invalid."
463  CALL flagerror(local_error,err,error,*999)
464  END SELECT
466  !u=exp(-kt)*sin(sqrt(k)*(x*cos(phi)+y*sin(phi)))
467  SELECT CASE(variable_type)
468  CASE(field_u_variable_type)
469  SELECT CASE(global_derivative)
470  CASE(no_global_deriv)
471  VALUE=exp(-k*time)*sin((sqrt(k))*(x(1)*cos(phi)+x(2)*sin(phi)))!Need to specify time, k and phi!
472  CASE(global_deriv_s1)
473  CALL flagerror("Not implemented.",err,error,*999)
474  CASE(global_deriv_s2)
475  CALL flagerror("Not implemented.",err,error,*999)
476  CASE(global_deriv_s1_s2)
477  CALL flagerror("Not implmented.",err,error,*999)
478  CASE DEFAULT
479  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
480  & " is invalid."
481  CALL flagerror(local_error,err,error,*999)
482  END SELECT
483  CASE(field_deludeln_variable_type)
484  SELECT CASE(global_derivative)
485  CASE(no_global_deriv)
486  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
487  CASE(global_deriv_s1)
488  CALL flagerror("Not implemented.",err,error,*999)
489  CASE(global_deriv_s2)
490  CALL flagerror("Not implemented.",err,error,*999)
491  CASE(global_deriv_s1_s2)
492  CALL flagerror("Not implemented.",err,error,*999)
493  CASE DEFAULT
494  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
495  & " is invalid."
496  CALL flagerror(local_error,err,error,*999)
497  END SELECT
498  CASE DEFAULT
499  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
500  & " is invalid."
501  CALL flagerror(local_error,err,error,*999)
502  END SELECT
504  !u=A1*exp(-t)*(x^2+y^2+z^2)
505  SELECT CASE(variable_type)
506  CASE(field_u_variable_type)
507  SELECT CASE(global_derivative)
508  CASE(no_global_deriv)
509  VALUE=a1*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
510  CASE(global_deriv_s1)
511  CALL flagerror("Not implemented.",err,error,*999)
512  CASE(global_deriv_s2)
513  CALL flagerror("Not implemented.",err,error,*999)
514  CASE(global_deriv_s1_s2)
515  CALL flagerror("Not implmented.",err,error,*999)
516  CASE DEFAULT
517  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
518  & " is invalid."
519  CALL flagerror(local_error,err,error,*999)
520  END SELECT
521  CASE(field_deludeln_variable_type)
522  SELECT CASE(global_derivative)
523  CASE(no_global_deriv)
524  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
525  CASE(global_deriv_s1)
526  CALL flagerror("Not implemented.",err,error,*999)
527  CASE(global_deriv_s2)
528  CALL flagerror("Not implemented.",err,error,*999)
529  CASE(global_deriv_s1_s2)
530  CALL flagerror("Not implemented.",err,error,*999)
531  CASE DEFAULT
532  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
533  & " is invalid."
534  CALL flagerror(local_error,err,error,*999)
535  END SELECT
536  CASE DEFAULT
537  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
538  & " is invalid."
539  CALL flagerror(local_error,err,error,*999)
540  END SELECT
541  CASE DEFAULT
542  local_error="The analytic function type of "// &
543  & trim(number_to_vstring(analytic_function_type,"*",err,error))// &
544  & " is invalid."
545  CALL flagerror(local_error,err,error,*999)
546  END SELECT
548  CALL flagerror("Not implemented.",err,error,*999)
550  SELECT CASE(analytic_function_type)
552  !u=exp(-kt)*sin(sqrt(k)*(x*cos(phi)+y*sin(phi)))
553  !These are parameters for the 3D analytical solution with a linear source
554  a = -0.25_dp
555  b = 0.5_dp
556  c = 0.5_dp
557  d = 0.5_dp
558  SELECT CASE(variable_type)
559  CASE(field_u_variable_type)
560  SELECT CASE(global_derivative)
561  CASE(no_global_deriv)
562  VALUE=exp(a*time)*exp(b*x(1))*exp(c*x(2))*exp(d*x(3))!Need to specify time, k and phi!
563  CASE(global_deriv_s1)
564  CALL flagerror("Not implemented.",err,error,*999)
565  CASE(global_deriv_s2)
566  CALL flagerror("Not implemented.",err,error,*999)
567  CASE(global_deriv_s1_s2)
568  CALL flagerror("Not implmented.",err,error,*999)
569  CASE DEFAULT
570  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
571  & " is invalid."
572  CALL flagerror(local_error,err,error,*999)
573  END SELECT
574  CASE(field_deludeln_variable_type)
575  SELECT CASE(global_derivative)
576  CASE(no_global_deriv)
577  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
578  CASE(global_deriv_s1)
579  CALL flagerror("Not implemented.",err,error,*999)
580  CASE(global_deriv_s2)
581  CALL flagerror("Not implemented.",err,error,*999)
582  CASE(global_deriv_s1_s2)
583  CALL flagerror("Not implemented.",err,error,*999)
584  CASE DEFAULT
585  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
586  & " is invalid."
587  CALL flagerror(local_error,err,error,*999)
588  END SELECT
589  CASE DEFAULT
590  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
591  & " is invalid."
592  CALL flagerror(local_error,err,error,*999)
593  END SELECT
594  CASE DEFAULT
595  local_error="The analytic function type of "// &
596  & trim(number_to_vstring(analytic_function_type,"*",err,error))// &
597  & " is invalid."
598  CALL flagerror(local_error,err,error,*999)
599  END SELECT
601  SELECT CASE(analytic_function_type)
603  !For del u/del t = del^2 u/del x^2 + a + bu + cu^m with a = 0 then
604  !u(x,t) = [+/-\beta + C.exp(\lamba.t+/\mu.x)]^(2/1-m) where
605  !\beta=\sqrt(-c/b); \lamba=b(1-m)(m+3)/(2(m+1)); \mu = \sqrt((b(1-m)^2)/(2.(m+1))
606  !see http://eqworld.ipmnet.ru/en/solutions/npde/npde1104.pdf
607  a_param=materials_parameters(1)
608  b_param=materials_parameters(2)
609  c_param=materials_parameters(3)
610  beta_param=sqrt(-c_param/b_param)
611  lambda_param=-5.0_dp*b_param/6.0_dp
612  mu_param=sqrt(b_param/6.0_dp)
613  const_param=1.0_dp
614  SELECT CASE(variable_type)
615  CASE(field_u_variable_type)
616  SELECT CASE(global_derivative)
617  CASE(no_global_deriv)
618  VALUE=1.0_dp/(beta_param+const_param*exp(lambda_param*time+mu_param*x(1)))**2
619  CASE(global_deriv_s1)
620  CALL flagerror("Not implemented.",err,error,*999)
621  CASE DEFAULT
622  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
623  & " is invalid."
624  CALL flagerror(local_error,err,error,*999)
625  END SELECT
626  CASE(field_deludeln_variable_type)
627  SELECT CASE(global_derivative)
628  CASE(no_global_deriv)
629  VALUE=0.0_dp
630  CASE(global_deriv_s1)
631  CALL flagerror("Not implemented.",err,error,*999)
632  CASE DEFAULT
633  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
634  & " is invalid."
635  CALL flagerror(local_error,err,error,*999)
636  END SELECT
637  CASE DEFAULT
638  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
639  & " is invalid."
640  CALL flagerror(local_error,err,error,*999)
641  END SELECT
642  CASE DEFAULT
643  local_error="The analytic function type of "// &
644  & trim(number_to_vstring(analytic_function_type,"*",err,error))// &
645  & " is invalid."
646  CALL flagerror(local_error,err,error,*999)
647  END SELECT
649  SELECT CASE(analytic_function_type)
651  !For del u/del t = del^2 u/del x^2 + a + b.e^c.u
652  !u(x,t) = -2/c.ln[+/-\beta + C.exp(+/-\mu.x-a.c.t/2)] where
653  !\beta=\sqrt(-b/a); \mu = \sqrt(a.c/2)
654  !see http://eqworld.ipmnet.ru/en/solutions/npde/npde1105.pdf
655  a_param=materials_parameters(1)
656  b_param=materials_parameters(2)
657  c_param=materials_parameters(3)
658  const_param=1.0_dp
659  beta_param=sqrt(-b_param/a_param)
660  mu_param=sqrt(a_param*c_param/2.0_dp)
661  SELECT CASE(variable_type)
662  CASE(field_u_variable_type)
663  SELECT CASE(global_derivative)
664  CASE(no_global_deriv)
665  VALUE=-2.0_dp/c_param*log(beta_param+const_param*exp(mu_param*x(1)-a_param*c_param*time/2.0_dp))
666  CASE(global_deriv_s1)
667  CALL flagerror("Not implemented.",err,error,*999)
668  CASE DEFAULT
669  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
670  & " is invalid."
671  CALL flagerror(local_error,err,error,*999)
672  END SELECT
673  CASE(field_deludeln_variable_type)
674  SELECT CASE(global_derivative)
675  CASE(no_global_deriv)
676  VALUE=0.0_dp
677  CASE(global_deriv_s1)
678  CALL flagerror("Not implemented.",err,error,*999)
679  CASE DEFAULT
680  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
681  & " is invalid."
682  CALL flagerror(local_error,err,error,*999)
683  END SELECT
684  CASE DEFAULT
685  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
686  & " is invalid."
687  CALL flagerror(local_error,err,error,*999)
688  END SELECT
689  CASE DEFAULT
690  local_error="The analytic function type of "// &
691  & trim(number_to_vstring(analytic_function_type,"*",err,error))// &
692  & " is invalid."
693  CALL flagerror(local_error,err,error,*999)
694  END SELECT
696  CALL flagerror("Not implemented.",err,error,*999)
698  CALL flagerror("Not implemented.",err,error,*999)
700  CALL flagerror("Not implemented.",err,error,*999)
702  CALL flagerror("Not implemented.",err,error,*999)
704  SELECT CASE(analytic_function_type)
706  SELECT CASE(variable_type)
707  CASE(field_u_variable_type)
708  SELECT CASE(global_derivative)
709  CASE(no_global_deriv)
710  VALUE=a1*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2))
711  CASE(global_deriv_s1)
712  CALL flagerror("Not implemented.",err,error,*999)
713  CASE(global_deriv_s2)
714  CALL flagerror("Not implemented.",err,error,*999)
715  CASE(global_deriv_s1_s2)
716  CALL flagerror("Not implmented.",err,error,*999)
717  CASE DEFAULT
718  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
719  & " is invalid."
720  CALL flagerror(local_error,err,error,*999)
721  END SELECT
722  CASE(field_deludeln_variable_type)
723  SELECT CASE(global_derivative)
724  CASE(no_global_deriv)
725  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
726  CASE(global_deriv_s1)
727  CALL flagerror("Not implemented.",err,error,*999)
728  CASE(global_deriv_s2)
729  CALL flagerror("Not implemented.",err,error,*999)
730  CASE(global_deriv_s1_s2)
731  CALL flagerror("Not implemented.",err,error,*999)
732  CASE DEFAULT
733  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
734  & " is invalid."
735  CALL flagerror(local_error,err,error,*999)
736  END SELECT
737  CASE(field_v_variable_type)
738  SELECT CASE(global_derivative)
739  CASE(no_global_deriv)
740  VALUE=a2*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2))
741  CASE(global_deriv_s1)
742  CALL flagerror("Not implemented.",err,error,*999)
743  CASE(global_deriv_s2)
744  CALL flagerror("Not implemented.",err,error,*999)
745  CASE(global_deriv_s1_s2)
746  CALL flagerror("Not implmented.",err,error,*999)
747  CASE DEFAULT
748  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
749  & " is invalid."
750  CALL flagerror(local_error,err,error,*999)
751  END SELECT
752  CASE(field_delvdeln_variable_type)
753  SELECT CASE(global_derivative)
754  CASE(no_global_deriv)
755  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
756  CASE(global_deriv_s1)
757  CALL flagerror("Not implemented.",err,error,*999)
758  CASE(global_deriv_s2)
759  CALL flagerror("Not implemented.",err,error,*999)
760  CASE(global_deriv_s1_s2)
761  CALL flagerror("Not implemented.",err,error,*999)
762  CASE DEFAULT
763  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
764  & " is invalid."
765  CALL flagerror(local_error,err,error,*999)
766  END SELECT
767  CASE DEFAULT
768  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
769  & " is invalid."
770  CALL flagerror(local_error,err,error,*999)
771  END SELECT
773  SELECT CASE(variable_type)
774  CASE(field_u_variable_type)
775  SELECT CASE(global_derivative)
776  CASE(no_global_deriv)
777  VALUE=a1*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
778  CASE(global_deriv_s1)
779  CALL flagerror("Not implemented.",err,error,*999)
780  CASE(global_deriv_s2)
781  CALL flagerror("Not implemented.",err,error,*999)
782  CASE(global_deriv_s1_s2)
783  CALL flagerror("Not implmented.",err,error,*999)
784  CASE DEFAULT
785  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
786  & " is invalid."
787  CALL flagerror(local_error,err,error,*999)
788  END SELECT
789  CASE(field_deludeln_variable_type)
790  SELECT CASE(global_derivative)
791  CASE(no_global_deriv)
792  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
793  CASE(global_deriv_s1)
794  CALL flagerror("Not implemented.",err,error,*999)
795  CASE(global_deriv_s2)
796  CALL flagerror("Not implemented.",err,error,*999)
797  CASE(global_deriv_s1_s2)
798  CALL flagerror("Not implemented.",err,error,*999)
799  CASE DEFAULT
800  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
801  & " is invalid."
802  CALL flagerror(local_error,err,error,*999)
803  END SELECT
804  CASE(field_v_variable_type)
805  SELECT CASE(global_derivative)
806  CASE(no_global_deriv)
807  VALUE=a2*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
808  CASE(global_deriv_s1)
809  CALL flagerror("Not implemented.",err,error,*999)
810  CASE(global_deriv_s2)
811  CALL flagerror("Not implemented.",err,error,*999)
812  CASE(global_deriv_s1_s2)
813  CALL flagerror("Not implmented.",err,error,*999)
814  CASE DEFAULT
815  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
816  & " is invalid."
817  CALL flagerror(local_error,err,error,*999)
818  END SELECT
819  CASE(field_delvdeln_variable_type)
820  SELECT CASE(global_derivative)
821  CASE(no_global_deriv)
822  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
823  CASE(global_deriv_s1)
824  CALL flagerror("Not implemented.",err,error,*999)
825  CASE(global_deriv_s2)
826  CALL flagerror("Not implemented.",err,error,*999)
827  CASE(global_deriv_s1_s2)
828  CALL flagerror("Not implemented.",err,error,*999)
829  CASE DEFAULT
830  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
831  & " is invalid."
832  CALL flagerror(local_error,err,error,*999)
833  END SELECT
834  CASE DEFAULT
835  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
836  & " is invalid."
837  CALL flagerror(local_error,err,error,*999)
838  END SELECT
840  SELECT CASE(variable_type)
841  CASE(field_u_variable_type)
842  SELECT CASE(global_derivative)
843  CASE(no_global_deriv)
844  VALUE=a1*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
845  CASE(global_deriv_s1)
846  CALL flagerror("Not implemented.",err,error,*999)
847  CASE(global_deriv_s2)
848  CALL flagerror("Not implemented.",err,error,*999)
849  CASE(global_deriv_s1_s2)
850  CALL flagerror("Not implmented.",err,error,*999)
851  CASE DEFAULT
852  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
853  & " is invalid."
854  CALL flagerror(local_error,err,error,*999)
855  END SELECT
856  CASE(field_deludeln_variable_type)
857  SELECT CASE(global_derivative)
858  CASE(no_global_deriv)
859  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
860  CASE(global_deriv_s1)
861  CALL flagerror("Not implemented.",err,error,*999)
862  CASE(global_deriv_s2)
863  CALL flagerror("Not implemented.",err,error,*999)
864  CASE(global_deriv_s1_s2)
865  CALL flagerror("Not implemented.",err,error,*999)
866  CASE DEFAULT
867  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
868  & " is invalid."
869  CALL flagerror(local_error,err,error,*999)
870  END SELECT
871  CASE(field_v_variable_type)
872  SELECT CASE(global_derivative)
873  CASE(no_global_deriv)
874  VALUE=a2*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
875  CASE(global_deriv_s1)
876  CALL flagerror("Not implemented.",err,error,*999)
877  CASE(global_deriv_s2)
878  CALL flagerror("Not implemented.",err,error,*999)
879  CASE(global_deriv_s1_s2)
880  CALL flagerror("Not implmented.",err,error,*999)
881  CASE DEFAULT
882  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
883  & " is invalid."
884  CALL flagerror(local_error,err,error,*999)
885  END SELECT
886  CASE(field_delvdeln_variable_type)
887  SELECT CASE(global_derivative)
888  CASE(no_global_deriv)
889  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
890  CASE(global_deriv_s1)
891  CALL flagerror("Not implemented.",err,error,*999)
892  CASE(global_deriv_s2)
893  CALL flagerror("Not implemented.",err,error,*999)
894  CASE(global_deriv_s1_s2)
895  CALL flagerror("Not implemented.",err,error,*999)
896  CASE DEFAULT
897  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
898  & " is invalid."
899  CALL flagerror(local_error,err,error,*999)
900  END SELECT
901  CASE(field_u1_variable_type)
902  SELECT CASE(global_derivative)
903  CASE(no_global_deriv)
904  VALUE=a3*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
905  CASE(global_deriv_s1)
906  CALL flagerror("Not implemented.",err,error,*999)
907  CASE(global_deriv_s2)
908  CALL flagerror("Not implemented.",err,error,*999)
909  CASE(global_deriv_s1_s2)
910  CALL flagerror("Not implmented.",err,error,*999)
911  CASE DEFAULT
912  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
913  & " is invalid."
914  CALL flagerror(local_error,err,error,*999)
915  END SELECT
916  CASE(field_delu1deln_variable_type)
917  SELECT CASE(global_derivative)
918  CASE(no_global_deriv)
919  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
920  CASE(global_deriv_s1)
921  CALL flagerror("Not implemented.",err,error,*999)
922  CASE(global_deriv_s2)
923  CALL flagerror("Not implemented.",err,error,*999)
924  CASE(global_deriv_s1_s2)
925  CALL flagerror("Not implemented.",err,error,*999)
926  CASE DEFAULT
927  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
928  & " is invalid."
929  CALL flagerror(local_error,err,error,*999)
930  END SELECT
931  CASE(field_u2_variable_type)
932  SELECT CASE(global_derivative)
933  CASE(no_global_deriv)
934  VALUE=a4*exp(-1.0_dp*time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
935  CASE(global_deriv_s1)
936  CALL flagerror("Not implemented.",err,error,*999)
937  CASE(global_deriv_s2)
938  CALL flagerror("Not implemented.",err,error,*999)
939  CASE(global_deriv_s1_s2)
940  CALL flagerror("Not implmented.",err,error,*999)
941  CASE DEFAULT
942  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
943  & " is invalid."
944  CALL flagerror(local_error,err,error,*999)
945  END SELECT
946  CASE(field_delu2deln_variable_type)
947  SELECT CASE(global_derivative)
948  CASE(no_global_deriv)
949  VALUE=0.0_dp!set to zero currently- actual value for diffusion solution needs adding
950  CASE(global_deriv_s1)
951  CALL flagerror("Not implemented.",err,error,*999)
952  CASE(global_deriv_s2)
953  CALL flagerror("Not implemented.",err,error,*999)
954  CASE(global_deriv_s1_s2)
955  CALL flagerror("Not implemented.",err,error,*999)
956  CASE DEFAULT
957  local_error="The global derivative index of "//trim(number_to_vstring(global_derivative,"*",err,error))// &
958  & " is invalid."
959  CALL flagerror(local_error,err,error,*999)
960  END SELECT
961  CASE DEFAULT
962  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
963  & " is invalid."
964  CALL flagerror(local_error,err,error,*999)
965  END SELECT
966  CASE DEFAULT
967  local_error="The analytic function type of "// &
968  & trim(number_to_vstring(analytic_function_type,"*",err,error))// &
969  & " is invalid."
970  CALL flagerror(local_error,err,error,*999)
971  END SELECT
972  CASE DEFAULT
973  local_error="The equations set subtype of "//trim(number_to_vstring(equations_subtype,"*",err,error))// &
974  & " is invalid."
975  CALL flagerror(local_error,err,error,*999)
976  END SELECT
977 
978  exits("Diffusion_AnalyticFunctionsEvaluate")
979  RETURN
980 999 errorsexits("Diffusion_AnalyticFunctionsEvaluate",err,error)
981  RETURN 1
983 
984  !
985  !================================================================================================================================
986  !
987 
989  SUBROUTINE diffusion_equation_equations_set_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
991  !Argument variables
992  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
993  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
994  INTEGER(INTG), INTENT(OUT) :: ERR
995  TYPE(varying_string), INTENT(OUT) :: ERROR
996  !Local Variables
997  TYPE(varying_string) :: LOCAL_ERROR
998 
999  enters("DIFFUSION_EQUATION_EQUATIONS_SET_SETUP",err,error,*999)
1000 
1001  IF(ASSOCIATED(equations_set)) THEN
1002  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
1003  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
1004  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
1005  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
1006  & err,error,*999)
1007  END IF
1008  SELECT CASE(equations_set%SPECIFICATION(3))
1010  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1012  !Need to define the functions diffusion_equation_equations_set_linear_source_setup etc
1013  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1015  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1017  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1019  !Need to define the functions diffusion_equation_equations_set_linear_source_setup etc
1020  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1022  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1024  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1026  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1028  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1030  CALL diffusion_equationssetlinearsetup(equations_set,equations_set_setup,err,error,*999)
1032  CALL diffusion_equationssetnonlinearsetup(equations_set,equations_set_setup,err,error,*999)
1034  CALL diffusion_equationssetnonlinearsetup(equations_set,equations_set_setup,err,error,*999)
1036  CALL flagerror("Not implemented.",err,error,*999)
1037 ! CALL Diffusion_EquationsSetNonlinearSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*999)
1039  CALL flagerror("Not implemented.",err,error,*999)
1040 ! CALL Diffusion_EquationsSetNonlinearSetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*999)
1041  CASE DEFAULT
1042  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1043  & " is not valid for a diffusion equation type of a classical field equation set class."
1044  CALL flagerror(local_error,err,error,*999)
1045  END SELECT
1046  ELSE
1047  CALL flagerror("Equations set is not associated.",err,error,*999)
1048  ENDIF
1049 
1050  exits("DIFFUSION_EQUATION_EQUATIONS_SET_SETUP")
1051  RETURN
1052 999 errorsexits("DIFFUSION_EQUATION_EQUATIONS_SET_SETUP",err,error)
1053  RETURN 1
1055 
1056  !
1057  !================================================================================================================================
1058  !
1059 
1061  SUBROUTINE diffusion_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)
1063  !Argument variables
1064  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1065  INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
1066  INTEGER(INTG), INTENT(OUT) :: ERR
1067  TYPE(varying_string), INTENT(OUT) :: ERROR
1068  !Local Variables
1069  TYPE(varying_string) :: LOCAL_ERROR
1070 
1071  enters("Diffusion_EquationsSetSolutionMethodSet",err,error,*999)
1072 
1073  IF(ASSOCIATED(equations_set)) THEN
1074  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
1075  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
1076  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
1077  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
1078  & err,error,*999)
1079  END IF
1080  SELECT CASE(equations_set%SPECIFICATION(3))
1089  SELECT CASE(solution_method)
1091  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
1093  CALL flagerror("Not implemented.",err,error,*999)
1095  CALL flagerror("Not implemented.",err,error,*999)
1097  CALL flagerror("Not implemented.",err,error,*999)
1099  CALL flagerror("Not implemented.",err,error,*999)
1101  CALL flagerror("Not implemented.",err,error,*999)
1102  CASE DEFAULT
1103  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
1104  CALL flagerror(local_error,err,error,*999)
1105  END SELECT
1106  CASE DEFAULT
1107  local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1108  & " is not valid for a diffusion equation type of an classical field equations set class."
1109  CALL flagerror(local_error,err,error,*999)
1110  END SELECT
1111  ELSE
1112  CALL flagerror("Equations set is not associated.",err,error,*999)
1113  ENDIF
1114 
1115  exits("Diffusion_EquationsSetSolutionMethodSet")
1116  RETURN
1117 999 errors("Diffusion_EquationsSetSolutionMethodSet",err,error)
1118  exits("Diffusion_EquationsSetSolutionMethodSet")
1119  RETURN 1
1120 
1122 
1123  !
1124  !================================================================================================================================
1125  !
1126 
1128  SUBROUTINE diffusion_equationssetspecificationset(equationsSet,specification,err,error,*)
1130  !Argument variables
1131  TYPE(equations_set_type), POINTER :: equationsSet
1132  INTEGER(INTG), INTENT(IN) :: specification(:)
1133  INTEGER(INTG), INTENT(OUT) :: err
1134  TYPE(varying_string), INTENT(OUT) :: error
1135  !Local Variables
1136  TYPE(varying_string) :: localError
1137  INTEGER(INTG) :: subtype
1138 
1139  enters("Diffusion_EquationsSetSpecificationSet",err,error,*999)
1140 
1141  IF(ASSOCIATED(equationsset)) THEN
1142  IF(SIZE(specification,1)/=3) THEN
1143  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
1144  & err,error,*999)
1145  END IF
1146  subtype=specification(3)
1147  SELECT CASE(subtype)
1162  !ok
1163  CASE DEFAULT
1164  localerror="The third equations set specification of "//trim(numbertovstring(subtype,"*",err,error))// &
1165  & " is not valid for a diffusion type of a classical field equations set."
1166  CALL flagerror(localerror,err,error,*999)
1167  END SELECT
1168  !Set full specification
1169  IF(ALLOCATED(equationsset%specification)) THEN
1170  CALL flagerror("Equations set specification is already allocated.",err,error,*999)
1171  ELSE
1172  ALLOCATE(equationsset%specification(3),stat=err)
1173  IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
1174  END IF
1175  equationsset%specification(1:3)=[equations_set_classical_field_class,equations_set_diffusion_equation_type,subtype]
1176  ELSE
1177  CALL flagerror("Equations set is not associated.",err,error,*999)
1178  END IF
1179 
1180  exits("Diffusion_EquationsSetSpecificationSet")
1181  RETURN
1182 999 errors("Diffusion_EquationsSetSpecificationSet",err,error)
1183  exits("Diffusion_EquationsSetSpecificationSet")
1184  RETURN 1
1185 
1187 
1188  !
1189  !================================================================================================================================
1190  !
1191 
1193  SUBROUTINE diffusion_equationssetlinearsetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
1195  !Argument variables
1196  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1197  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
1198  INTEGER(INTG), INTENT(OUT) :: ERR
1199  TYPE(varying_string), INTENT(OUT) :: ERROR
1200  !Local Variables
1201  INTEGER(INTG) :: component_idx,GEOMETRIC_MESH_COMPONENT,GEOMETRIC_SCALING_TYPE,NUMBER_OF_ANALYTIC_COMPONENTS, &
1202  & NUMBER_OF_DIMENSIONS, NUMBER_OF_MATERIALS_COMPONENTS, NUMBER_OF_SOURCE_COMPONENTS,imy_matrix,Ncompartments, &
1203  & GEOMETRIC_COMPONENT_NUMBER
1204  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
1205  TYPE(equations_type), POINTER :: EQUATIONS
1206  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
1207  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
1208  TYPE(equations_set_analytic_type), POINTER :: EQUATIONS_ANALYTIC
1209  TYPE(equations_set_equations_set_field_type), POINTER :: EQUATIONS_EQUATIONS_SET_FIELD
1210  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
1211  TYPE(equations_set_source_type), POINTER :: EQUATIONS_SOURCE
1212  TYPE(field_type), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,GEOMETRIC_FIELD,EQUATIONS_SET_FIELD_FIELD
1213  TYPE(varying_string) :: LOCAL_ERROR
1214  INTEGER(INTG) :: num_var,num_var_count,NUMBER_OF_MATERIALS_COUPLING_COMPONENTS
1215  INTEGER(INTG) :: EQUATIONS_SET_FIELD_NUMBER_OF_VARIABLES,EQUATIONS_SET_FIELD_NUMBER_OF_COMPONENTS
1216  INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
1217  INTEGER(INTG), ALLOCATABLE :: VARIABLE_TYPES(:),VARIABLE_U_TYPES(:),COUPLING_MATRIX_STORAGE_TYPE(:), &
1218  & COUPLING_MATRIX_STRUCTURE_TYPE(:)
1219 
1220  enters("Diffusion_EquationsSetLinearSetup",err,error,*999)
1221 
1222  NULLIFY(equations)
1223  NULLIFY(equations_mapping)
1224  NULLIFY(equations_matrices)
1225  NULLIFY(geometric_decomposition)
1226 
1227  IF(ASSOCIATED(equations_set)) THEN
1228  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
1229  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
1230  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
1231  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
1232  & err,error,*999)
1233  END IF
1234  IF(equations_set%SPECIFICATION(3)==equations_set_no_source_diffusion_subtype .OR. &
1235  & equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
1236  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
1237  & equations_set%SPECIFICATION(3)==equations_set_no_source_ale_diffusion_subtype .OR. &
1238  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
1239  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
1240  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype .OR. &
1241  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype .OR. &
1242  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_advec_diff_subtype .OR. &
1243  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_advec_diff_supg_subtype) THEN
1244  SELECT CASE(equations_set_setup%SETUP_TYPE)
1246  SELECT CASE(equations_set_setup%ACTION_TYPE)
1249  & err,error,*999)
1250  IF(equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype .OR. &
1251  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_advec_diff_subtype .OR. &
1252  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_advec_diff_supg_subtype) THEN
1253  equations_set_field_number_of_variables = 1
1254  equations_set_field_number_of_components = 2
1255  equations_equations_set_field=>equations_set%EQUATIONS_SET_FIELD
1256  IF(equations_equations_set_field%EQUATIONS_SET_FIELD_AUTO_CREATED) THEN
1257  !Create the auto created equations set field
1258  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
1259  & equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,err,error,*999)
1260  CALL field_label_set(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,"Equations Set Field",err,error,*999)
1261  CALL field_type_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,field_general_type,&
1262  & err,error,*999)
1263  CALL field_dependent_type_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,&
1264  & field_independent_type,err,error,*999)
1265  CALL field_number_of_variables_set(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD, &
1266  & equations_set_field_number_of_variables,err,error,*999)
1267  CALL field_variable_types_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,&
1268  & [field_u_variable_type],err,error,*999)
1269  CALL field_variable_label_set(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,field_u_variable_type, &
1270  & "Equations",err,error,*999)
1271  CALL field_dimension_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,field_u_variable_type, &
1272  & field_vector_dimension_type,err,error,*999)
1273  CALL field_data_type_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,field_u_variable_type, &
1274  & field_intg_type,err,error,*999)
1275  CALL field_number_of_components_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,&
1276  & field_u_variable_type,equations_set_field_number_of_components,err,error,*999)
1277  ELSE
1278  !Check the user specified field
1279  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1280  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1281  CALL field_number_of_variables_check(equations_set_setup%FIELD,equations_set_field_number_of_variables, &
1282  & err,error,*999)
1283  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
1284  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1285  & err,error,*999)
1286  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_intg_type,err,error,*999)
1287  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
1288  & equations_set_field_number_of_components,err,error,*999)
1289  ENDIF
1290  ENDIF
1292  IF(equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype .OR. &
1293  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_advec_diff_subtype .OR. &
1294  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_advec_diff_supg_subtype) THEN
1295  IF(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_AUTO_CREATED) THEN
1296  CALL field_create_finish(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,err,error,*999)
1297  CALL field_component_values_initialise(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,&
1298  & field_u_variable_type,field_values_set_type, 1, 1_intg, err, error, *999)
1299  CALL field_component_values_initialise(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,&
1300  & field_u_variable_type,field_values_set_type, 2, 1_intg, err, error, *999)
1301  ENDIF
1302  ENDIF
1303 !!TODO: Check valid setup
1304  CASE DEFAULT
1305  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1306  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1307  & " is invalid for a linear diffusion equation."
1308  CALL flagerror(local_error,err,error,*999)
1309  END SELECT
1311  SELECT CASE(equations_set%SPECIFICATION(3))
1316  !do nothing
1318  SELECT CASE(equations_set_setup%ACTION_TYPE)
1320  equations_set_field_number_of_components = 2
1321  equations_equations_set_field=>equations_set%EQUATIONS_SET_FIELD
1322  IF(equations_equations_set_field%EQUATIONS_SET_FIELD_AUTO_CREATED) THEN
1323  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1324  CALL field_mesh_decomposition_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,&
1325  & geometric_decomposition,err,error,*999)
1326  CALL field_geometric_field_set_and_lock(equations_equations_set_field%EQUATIONS_SET_FIELD_FIELD,&
1327  & equations_set%GEOMETRY%GEOMETRIC_FIELD,err,error,*999)
1328  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1329  & 1,geometric_component_number,err,error,*999)
1330  DO component_idx = 1, equations_set_field_number_of_components
1331  CALL field_component_mesh_component_set_and_lock(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD, &
1332  & field_u_variable_type,component_idx,geometric_component_number,err,error,*999)
1333  CALL field_component_interpolation_set_and_lock(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD, &
1334  & field_u_variable_type,component_idx,field_constant_interpolation,err,error,*999)
1335  END DO
1336  !Default the field scaling to that of the geometric field
1337  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1338  CALL field_scaling_type_set(equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD,geometric_scaling_type, &
1339  & err,error,*999)
1340  ELSE
1341  !Do nothing
1342  ENDIF
1344  !Do nothing
1345  CASE DEFAULT
1346  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1347  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1348  & " is invalid for a linear diffusion equation."
1349  CALL flagerror(local_error,err,error,*999)
1350  END SELECT
1354  SELECT CASE(equations_set_setup%ACTION_TYPE)
1356  CALL field_parameter_set_create(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
1357  & field_mesh_displacement_set_type, err, error, *999)
1358  CALL field_parameter_set_create(equations_set%GEOMETRY%GEOMETRIC_FIELD, field_u_variable_type, &
1359  & field_mesh_velocity_set_type, err, error, *999)
1361  !Do nothing
1362  CASE DEFAULT
1363  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1364  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1365  & " is invalid for a linear diffusion equation."
1366  CALL flagerror(local_error,err,error,*999)
1367  END SELECT
1368  END SELECT
1369  !-----------------------------------------------------------------
1370  ! D e p e n d e n t f i e l d
1371  !-----------------------------------------------------------------
1373  SELECT CASE(equations_set_setup%ACTION_TYPE)
1375  SELECT CASE(equations_set%SPECIFICATION(3))
1384  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
1385  !Create the auto created dependent field
1386  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
1387  & dependent_field,err,error,*999)
1388  CALL field_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,"Dependent Field",err,error,*999)
1389  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
1390  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
1391  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1392  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
1393  & err,error,*999)
1394  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
1395  & geometric_field,err,error,*999)
1396  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
1397  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,[field_u_variable_type, &
1398  & field_deludeln_variable_type],err,error,*999)
1399  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
1400  & "U",err,error,*999)
1401  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
1402  & "del U/del n",err,error,*999)
1403  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
1404  & field_scalar_dimension_type,err,error,*999)
1405  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
1406  & field_scalar_dimension_type,err,error,*999)
1407  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
1408  & field_dp_type,err,error,*999)
1409  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
1410  & field_dp_type,err,error,*999)
1411  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
1412  & err,error,*999)
1413  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
1414  & field_deludeln_variable_type,1,err,error,*999)
1415  !Default to the geometric interpolation setup
1416  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
1417  & geometric_mesh_component,err,error,*999)
1418  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
1419  & geometric_mesh_component,err,error,*999)
1420  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
1421  & geometric_mesh_component,err,error,*999)
1422  SELECT CASE(equations_set%SOLUTION_METHOD)
1424  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
1425  & field_u_variable_type,1,field_node_based_interpolation,err,error,*999)
1426  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
1427  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
1428  !Default the scaling to the geometric field scaling
1429  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1430  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
1432  CALL flagerror("Not implemented.",err,error,*999)
1434  CALL flagerror("Not implemented.",err,error,*999)
1436  CALL flagerror("Not implemented.",err,error,*999)
1438  CALL flagerror("Not implemented.",err,error,*999)
1440  CALL flagerror("Not implemented.",err,error,*999)
1441  CASE DEFAULT
1442  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1443  & " is invalid."
1444  CALL flagerror(local_error,err,error,*999)
1445  END SELECT
1446  ELSE
1447  SELECT CASE(equations_set%SPECIFICATION(3))
1449  !Check the field created by advection-diffusion routines for the coupled problem
1450  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1451  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
1452  CALL field_number_of_variables_check(equations_set_setup%FIELD,4,err,error,*999)
1453  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type,field_deludeln_variable_type, &
1454  & field_v_variable_type,field_delvdeln_variable_type],err,error,*999)
1455  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type, &
1456  & field_scalar_dimension_type,err,error,*999)
1457  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
1458  & err,error,*999)
1459  CALL field_dimension_check(equations_set_setup%FIELD,field_v_variable_type, &
1460  & field_scalar_dimension_type,err,error,*999)
1461  CALL field_dimension_check(equations_set_setup%FIELD,field_delvdeln_variable_type,field_scalar_dimension_type, &
1462  & err,error,*999)
1463  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1464  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
1465  CALL field_data_type_check(equations_set_setup%FIELD,field_v_variable_type,field_dp_type,err,error,*999)
1466  CALL field_data_type_check(equations_set_setup%FIELD,field_delvdeln_variable_type,field_dp_type,err,error,*999)
1467  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1468  & number_of_dimensions,err,error,*999)
1469  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1, &
1470  & err,error,*999)
1471  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
1472  & err,error,*999)
1473  CALL field_number_of_components_check(equations_set_setup%FIELD,field_v_variable_type,1, &
1474  & err,error,*999)
1475  CALL field_number_of_components_check(equations_set_setup%FIELD,field_delvdeln_variable_type,1, &
1476  & err,error,*999)
1477  SELECT CASE(equations_set%SOLUTION_METHOD)
1479  component_idx=1
1480  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,component_idx, &
1481  & field_node_based_interpolation,err,error,*999)
1482  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
1483  & component_idx,field_node_based_interpolation,err,error,*999)
1484  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_v_variable_type,component_idx, &
1485  & field_node_based_interpolation,err,error,*999)
1486  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_delvdeln_variable_type, &
1487  & component_idx,field_node_based_interpolation,err,error,*999)
1489  CALL flagerror("Not implemented.",err,error,*999)
1491  CALL flagerror("Not implemented.",err,error,*999)
1493  CALL flagerror("Not implemented.",err,error,*999)
1495  CALL flagerror("Not implemented.",err,error,*999)
1497  CALL flagerror("Not implemented.",err,error,*999)
1498  CASE DEFAULT
1499  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1500  & " is invalid."
1501  CALL flagerror(local_error,err,error,*999)
1502  END SELECT
1503 
1505  !uses number of compartments to check that appropriate number and type of variables have been set on the
1506  !dependent field
1507  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1508  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
1509  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1510  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1511  & field_values_set_type,equations_set_field_data,err,error,*999)
1512  ncompartments=equations_set_field_data(2)
1513  CALL field_number_of_variables_check(equations_set_setup%FIELD,2*ncompartments,err,error,*999)
1514  !Create & populate array storing all of the relevant variable types against which to check the field variables
1515  ALLOCATE(variable_types(2*ncompartments))
1516  DO num_var=1,ncompartments
1517  variable_types(2*num_var-1)=field_u_variable_type+(field_number_of_variable_subtypes*(num_var-1))
1518  variable_types(2*num_var)=field_deludeln_variable_type+(field_number_of_variable_subtypes*(num_var-1))
1519  ENDDO
1520  CALL field_variable_types_check(equations_set_setup%FIELD,variable_types,err,error,*999)
1521 
1522  DO num_var=1,2*ncompartments
1523  CALL field_dimension_check(equations_set_setup%FIELD,variable_types(num_var), &
1524  & field_scalar_dimension_type,err,error,*999)
1525  CALL field_data_type_check(equations_set_setup%FIELD,variable_types(num_var),field_dp_type,err,error,*999)
1526  CALL field_number_of_components_check(equations_set_setup%FIELD,variable_types(num_var),1, &
1527  & err,error,*999)
1528  ENDDO
1529  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1530  & number_of_dimensions,err,error,*999)
1531  SELECT CASE(equations_set%SOLUTION_METHOD)
1533  component_idx=1
1534  DO num_var=1,2*ncompartments
1535  CALL field_component_interpolation_check(equations_set_setup%FIELD,variable_types(num_var),component_idx, &
1536  & field_node_based_interpolation,err,error,*999)
1537  ENDDO
1539  CALL flagerror("Not implemented.",err,error,*999)
1541  CALL flagerror("Not implemented.",err,error,*999)
1543  CALL flagerror("Not implemented.",err,error,*999)
1545  CALL flagerror("Not implemented.",err,error,*999)
1547  CALL flagerror("Not implemented.",err,error,*999)
1548  CASE DEFAULT
1549  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1550  & " is invalid."
1551  CALL flagerror(local_error,err,error,*999)
1552  END SELECT
1555  CALL flagerror("Not implemented.",err,error,*999)
1556  CASE DEFAULT
1557  !Check the user specified field
1558  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1559  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
1560  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
1561  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type, &
1562  & field_deludeln_variable_type],err,error,*999)
1563  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type, &
1564  & err,error,*999)
1565  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
1566  & err,error,*999)
1567  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1568  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
1569  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1570  & number_of_dimensions,err,error,*999)
1571  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions, &
1572  & err,error,*999)
1573  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
1574  & number_of_dimensions,err,error,*999)
1575  SELECT CASE(equations_set%SOLUTION_METHOD)
1577  DO component_idx=1,number_of_dimensions
1578  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,component_idx, &
1579  & field_node_based_interpolation,err,error,*999)
1580  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
1581  & component_idx,field_node_based_interpolation,err,error,*999)
1582  ENDDO !component_idx
1584  CALL flagerror("Not implemented.",err,error,*999)
1586  CALL flagerror("Not implemented.",err,error,*999)
1588  CALL flagerror("Not implemented.",err,error,*999)
1590  CALL flagerror("Not implemented.",err,error,*999)
1592  CALL flagerror("Not implemented.",err,error,*999)
1593  CASE DEFAULT
1594  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1595  & " is invalid."
1596  CALL flagerror(local_error,err,error,*999)
1597  END SELECT
1598  END SELECT
1599  ENDIF
1600  END SELECT
1602  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
1603  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
1604  ENDIF
1605  CASE DEFAULT
1606  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1607  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1608  & " is invalid for a linear diffusion equation"
1609  CALL flagerror(local_error,err,error,*999)
1610  END SELECT
1611  !-----------------------------------------------------------------
1612  ! M a t e r i a l s f i e l d
1613  !-----------------------------------------------------------------
1615  SELECT CASE(equations_set_setup%ACTION_TYPE)
1617  equations_materials=>equations_set%MATERIALS
1618  IF(ASSOCIATED(equations_materials)) THEN
1619  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
1620  !Create the auto created materials field
1621  SELECT CASE(equations_set%SPECIFICATION(3))
1629  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
1630  & materials_field,err,error,*999)
1631  CALL field_label_set(equations_materials%MATERIALS_FIELD,"Materials Field",err,error,*999)
1632  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
1633  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
1634  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1635  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
1636  & err,error,*999)
1637  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
1638  & geometric_field,err,error,*999)
1639  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
1640  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type], &
1641  & err,error,*999)
1642  CALL field_variable_label_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1643  & "Materials",err,error,*999)
1644  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1645  & field_vector_dimension_type,err,error,*999)
1646  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1647  & field_dp_type,err,error,*999)
1648  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1649  & number_of_dimensions,err,error,*999)
1650  IF(equations_set%SPECIFICATION(3)==equations_set_no_source_diffusion_subtype .OR. &
1651  & equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
1652  & equations_set%SPECIFICATION(3)==equations_set_no_source_ale_diffusion_subtype .OR. &
1653  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype) THEN
1654  number_of_materials_components=number_of_dimensions
1655  ELSEIF(equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
1656  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype) THEN
1657  !Linear source. Materials field components are 1 for each dimension and 1 for the linear source
1658  !i.e., k and a in div(k.grad(u(x)))=a(x)u(x)+c(x)
1659  number_of_materials_components=number_of_dimensions+1
1660  ELSE
1661  number_of_materials_components=number_of_dimensions+2
1662  ENDIF
1663  !Set the number of materials components
1664  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1665  & number_of_materials_components,err,error,*999)
1666  !Default the k materials components to the geometric interpolation setup with constant interpolation
1667  DO component_idx=1,number_of_dimensions
1668  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1669  & component_idx,geometric_mesh_component,err,error,*999)
1670  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1671  & component_idx,field_constant_interpolation,err,error,*999)
1672  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1673  & component_idx,geometric_mesh_component,err,error,*999)
1674  ENDDO !component_idx
1675  !Default the source materials components to the first component geometric interpolation with constant
1676  !interpolation
1677  IF(equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
1678  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype) THEN
1679  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1680  & 1,geometric_mesh_component,err,error,*999)
1681  DO component_idx=number_of_dimensions+1,number_of_materials_components
1682  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1683  & component_idx,geometric_mesh_component,err,error,*999)
1684  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1685  & component_idx,field_constant_interpolation,err,error,*999)
1686  ENDDO !components_idx
1687  ENDIF
1688  !Default the field scaling to that of the geometric field
1689  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1690  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
1692  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
1693  & materials_field,err,error,*999)
1694  CALL field_label_set(equations_materials%MATERIALS_FIELD,"Materials Field",err,error,*999)
1695  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
1696  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
1697  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1698  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
1699  & err,error,*999)
1700  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
1701  & geometric_field,err,error,*999)
1702  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,2,err,error,*999)
1703  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type, &
1704  & field_v_variable_type],err,error,*999)
1705  CALL field_variable_label_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1706  & "Materials",err,error,*999)
1707  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1708  & field_vector_dimension_type,err,error,*999)
1709  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1710  & field_dp_type,err,error,*999)
1711  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1712  & field_vector_dimension_type,err,error,*999)
1713  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1714  & field_dp_type,err,error,*999)
1715  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1716  & number_of_dimensions,err,error,*999)
1717  number_of_materials_components=number_of_dimensions
1718  !Set the number of materials components
1719  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1720  & number_of_materials_components,err,error,*999)
1721  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1722  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1723  & field_values_set_type,equations_set_field_data,err,error,*999)
1724  ncompartments=equations_set_field_data(2)
1725  number_of_materials_coupling_components=ncompartments
1726  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1727  & number_of_materials_coupling_components,err,error,*999)
1728  !Default the k materials components to the geometric interpolation setup with constant interpolation
1729  DO component_idx=1,number_of_dimensions
1730  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1731  & component_idx,geometric_mesh_component,err,error,*999)
1732  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1733  & component_idx,field_constant_interpolation,err,error,*999)
1734  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1735  & component_idx,geometric_mesh_component,err,error,*999)
1736  ENDDO !component_idx
1737  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1738  & 1,geometric_mesh_component,err,error,*999)
1739  DO component_idx=1,number_of_materials_coupling_components
1740  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1741  & component_idx,field_constant_interpolation,err,error,*999)
1742  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1743  & component_idx,geometric_mesh_component,err,error,*999)
1744  ENDDO
1745  !Default the field scaling to that of the geometric field
1746  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1747  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
1748  END SELECT
1749  ELSE
1750  !Check the user specified field
1751  SELECT CASE(equations_set%SPECIFICATION(3))
1759 
1760  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
1761  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1762  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
1763  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
1764  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1765  & err,error,*999)
1766  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1767  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1768  & number_of_dimensions,err,error,*999)
1769  IF(equations_set%SPECIFICATION(3)==equations_set_no_source_diffusion_subtype .OR. &
1770  & equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
1771  & equations_set%SPECIFICATION(3)==equations_set_no_source_ale_diffusion_subtype .OR. &
1772  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype) THEN
1773  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions, &
1774  & err,error,*999)
1775  ELSE
1776  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+1, &
1777  & err,error,*999)
1778  ENDIF
1780  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
1781  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1782  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
1783  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type, &
1784  & field_v_variable_type],err,error,*999)
1785  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1786  & err,error,*999)
1787  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1788  CALL field_dimension_check(equations_set_setup%FIELD,field_v_variable_type,field_vector_dimension_type, &
1789  & err,error,*999)
1790  CALL field_data_type_check(equations_set_setup%FIELD,field_v_variable_type,field_dp_type,err,error,*999)
1791  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1792  & number_of_dimensions,err,error,*999)
1793  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions, &
1794  & err,error,*999)
1795  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1796  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1797  & field_values_set_type,equations_set_field_data,err,error,*999)
1798  ncompartments=equations_set_field_data(2)
1799  number_of_materials_coupling_components=ncompartments
1800  CALL field_number_of_components_check(equations_set_setup%FIELD,field_v_variable_type, &
1801  & number_of_materials_coupling_components,err,error,*999)
1802  END SELECT
1803  ENDIF
1804  ELSE
1805  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1806  ENDIF
1808  equations_materials=>equations_set%MATERIALS
1809  IF(ASSOCIATED(equations_materials)) THEN
1810  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
1811  !Finish creating the materials field
1812  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
1813  !Set the default values for the materials field
1814  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1815  & number_of_dimensions,err,error,*999)
1816  IF(equations_set%SPECIFICATION(3)==equations_set_no_source_diffusion_subtype .OR. &
1817  & equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
1818  & equations_set%SPECIFICATION(3)==equations_set_no_source_ale_diffusion_subtype .OR. &
1819  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype) THEN
1820  number_of_materials_components=number_of_dimensions
1821  ELSE IF(equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
1822  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype) THEN
1823  !Constant source. Materials field components are 1 for each dimension and 1 for the constant source
1824  !i.e., k and c in div(k.grad(u(x)))=c(x)
1825  number_of_materials_components=number_of_dimensions+1
1826  ENDIF
1827  !First set the k values to 1.0
1828  DO component_idx=1,number_of_dimensions
1829  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1830  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
1831  ENDDO !component_idx
1832  IF(equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype) THEN
1833  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
1834  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
1835  & field_values_set_type,equations_set_field_data,err,error,*999)
1836  ncompartments=equations_set_field_data(2)
1837  DO component_idx=1,ncompartments
1838  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_v_variable_type, &
1839  & field_values_set_type,component_idx,0.0_dp,err,error,*999)
1840  ENDDO !component_idx
1841  ENDIF
1842  IF(equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
1843  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype) THEN
1844  !Now set the linear source values to 1.0
1845  DO component_idx=number_of_dimensions+1,number_of_materials_components
1846  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1847  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
1848  ENDDO !component_idx
1849  ENDIF
1850  ENDIF
1851  ELSE
1852  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1853  ENDIF
1854  CASE DEFAULT
1855  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1856  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1857  & " is invalid for a linear diffusion equation."
1858  CALL flagerror(local_error,err,error,*999)
1859  END SELECT
1860  !-----------------------------------------------------------------
1861  ! S o u r c e f i e l d
1862  !-----------------------------------------------------------------
1864  SELECT CASE(equations_set_setup%ACTION_TYPE)
1866  equations_source=>equations_set%SOURCE
1867  IF(ASSOCIATED(equations_source)) THEN
1868  IF(equations_source%SOURCE_FIELD_AUTO_CREATED) THEN
1869  !Create the auto created source field
1870  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_source% &
1871  & source_field,err,error,*999)
1872  CALL field_label_set(equations_source%SOURCE_FIELD,"Source Field",err,error,*999)
1873  CALL field_type_set_and_lock(equations_source%SOURCE_FIELD,field_general_type,err,error,*999)
1874  CALL field_dependent_type_set_and_lock(equations_source%SOURCE_FIELD,field_independent_type,err,error,*999)
1875  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1876  CALL field_mesh_decomposition_set_and_lock(equations_source%SOURCE_FIELD,geometric_decomposition, &
1877  & err,error,*999)
1878  CALL field_geometric_field_set_and_lock(equations_source%SOURCE_FIELD,equations_set%GEOMETRY% &
1879  & geometric_field,err,error,*999)
1880  CALL field_number_of_variables_set_and_lock(equations_source%SOURCE_FIELD,1,err,error,*999)
1881  CALL field_variable_types_set_and_lock(equations_source%SOURCE_FIELD,[field_u_variable_type], &
1882  & err,error,*999)
1883  CALL field_variable_label_set(equations_source%SOURCE_FIELD,field_u_variable_type, &
1884  & "Source",err,error,*999)
1885  CALL field_dimension_set_and_lock(equations_source%SOURCE_FIELD,field_u_variable_type, &
1886  & field_scalar_dimension_type,err,error,*999)
1887  CALL field_data_type_set_and_lock(equations_source%SOURCE_FIELD,field_u_variable_type, &
1888  & field_dp_type,err,error,*999)
1889  number_of_source_components=1
1890  !Set the number of source components
1891  CALL field_number_of_components_set_and_lock(equations_source%SOURCE_FIELD,field_u_variable_type, &
1892  & number_of_source_components,err,error,*999)
1893  !Default the source components to the geometric interpolation setup with constant interpolation
1894  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
1895  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
1896  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
1897  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
1898  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype .OR. &
1899  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype) THEN
1900  DO component_idx=1,number_of_source_components
1901  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1902  & component_idx,geometric_mesh_component,err,error,*999)
1903  CALL field_component_mesh_component_set(equations_source%SOURCE_FIELD,field_u_variable_type, &
1904  & component_idx,geometric_mesh_component,err,error,*999)
1905  CALL field_component_interpolation_set(equations_source%SOURCE_FIELD,field_u_variable_type, &
1906  & component_idx,field_node_based_interpolation,err,error,*999)
1907  ENDDO !component_idx
1908  ENDIF
1909  !Default the field scaling to that of the geometric field
1910  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1911  CALL field_scaling_type_set(equations_source%SOURCE_FIELD,geometric_scaling_type,err,error,*999)
1912  ELSE
1913  !Check the user specified field
1914  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1915  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1916  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
1917  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
1918  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type, &
1919  & err,error,*999)
1920  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1921  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1, &
1922  & err,error,*999)
1923  ENDIF
1924  ELSE
1925  CALL flagerror("Equations set source is not associated.",err,error,*999)
1926  ENDIF
1928  equations_source=>equations_set%SOURCE
1929  IF(ASSOCIATED(equations_source)) THEN
1930  IF(equations_source%SOURCE_FIELD_AUTO_CREATED) THEN
1931  !Finish creating the source field
1932  CALL field_create_finish(equations_source%SOURCE_FIELD,err,error,*999)
1933  !Set the default values for the source field
1934  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1935  & number_of_dimensions,err,error,*999)
1936  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
1937  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
1938  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
1939  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
1940  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype .OR. &
1941  & equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype) THEN
1942  number_of_source_components=1
1943  ELSE
1944  number_of_source_components=0
1945  ENDIF
1946  !Now set the source values to 1.0
1947  DO component_idx=1,number_of_source_components
1948  CALL field_component_values_initialise(equations_source%SOURCE_FIELD,field_u_variable_type, &
1949  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
1950  ENDDO !component_idx
1951  ENDIF
1952  ELSE
1953  CALL flagerror("Equations set source is not associated.",err,error,*999)
1954  ENDIF
1955  CASE DEFAULT
1956  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1957  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1958  & " is invalid for a linear diffusion equation."
1959  CALL flagerror(local_error,err,error,*999)
1960  END SELECT
1961  !-----------------------------------------------------------------
1962  ! A n a l y t i c T y p e
1963  !-----------------------------------------------------------------
1965  SELECT CASE(equations_set_setup%ACTION_TYPE)
1967  equations_analytic=>equations_set%ANALYTIC
1968  IF(ASSOCIATED(equations_analytic)) THEN
1969  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
1970  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
1971  IF(ASSOCIATED(dependent_field)) THEN
1972  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
1973  IF(ASSOCIATED(geometric_field)) THEN
1974  equations_materials=>equations_set%MATERIALS
1975  IF(ASSOCIATED(equations_materials)) THEN
1976  IF(equations_materials%MATERIALS_FINISHED) THEN
1977  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions, &
1978  & err,error,*999)
1979  equations_set%ANALYTIC%ANALYTIC_USER_PARAMS(1)=0.0_dp
1980  SELECT CASE(equations_set%SPECIFICATION(3))
1982  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
1984  IF(number_of_dimensions/=1) THEN
1985  local_error="The number of geometric dimensions of "// &
1986  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
1987  & " is invalid. The analytic function type of "// &
1988  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1989  & " for a no source diffusion equation requires that there be 1 geometric dimension."
1990  CALL flagerror(local_error,err,error,*999)
1991  ENDIF
1992  !Check the materials values are constant
1993  CALL field_component_interpolation_check(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1994  & 1,field_constant_interpolation,err,error,*999)
1995  !Set number of analytic field components
1996  number_of_analytic_components=4
1997  !Set analytic function type
1998  equations_analytic%ANALYTIC_FUNCTION_TYPE=equations_set_diffusion_equation_one_dim_1
2000  !Check that domain is 2D
2001  IF(number_of_dimensions/=2) THEN
2002  local_error="The number of geometric dimensions of "// &
2003  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2004  & " is invalid. The analytic function type of "// &
2005  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2006  & " for a no source diffusion equation requires that there be 2 geometric dimensions."
2007  CALL flagerror(local_error,err,error,*999)
2008  ENDIF
2009  !Set number of analytic field components
2010  number_of_analytic_components=0
2011  !Set analytic function type
2012  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_diffusion_equation_two_dim_1
2013  CASE DEFAULT
2014  local_error="The specified analytic function type of "// &
2015  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2016  & " is invalid for a no source diffusion equation."
2017  CALL flagerror(local_error,err,error,*999)
2018  END SELECT
2020  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
2022  !Check that domain is 3D
2023  IF(number_of_dimensions/=3) THEN
2024  local_error="The number of geometric dimensions of "// &
2025  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2026  & " is invalid. The analytic function type of "// &
2027  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2028  & " for a constant source diffusion equation requires that there be 3 geometric dimensions."
2029  CALL flagerror(local_error,err,error,*999)
2030  ENDIF
2031  !Set number of analytic field components
2032  number_of_analytic_components=0
2033  !Set analytic function type
2034  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_diffusion_equation_three_dim_1
2035  CASE DEFAULT
2036  local_error="The specified analytic function type of "// &
2037  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2038  & " is invalid for a constant source diffusion equation."
2039  CALL flagerror(local_error,err,error,*999)
2040  END SELECT
2042  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
2044  !Check that domain is 3D
2045  IF(number_of_dimensions/=3) THEN
2046  local_error="The number of geometric dimensions of "// &
2047  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2048  & " is invalid. The analytic function type of "// &
2049  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2050  & " for a linear source diffusion equation requires that there be 3 geometric dimensions."
2051  CALL flagerror(local_error,err,error,*999)
2052  ENDIF
2053  !Set number of analytic field components
2054  number_of_analytic_components=0
2055  !Set analytic function type
2056  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE= &
2058  CASE DEFAULT
2059  local_error="The specified analytic function type of "// &
2060  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2061  & " is invalid for a linear source diffusion equation."
2062  CALL flagerror(local_error,err,error,*999)
2063  END SELECT
2065  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
2067  !Check that domain is 2D
2068  IF(number_of_dimensions/=2) THEN
2069  local_error="The number of geometric dimensions of "// &
2070  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2071  & " is invalid. The analytic function type of "// &
2072  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2073  & " for a multi-compartment diffusion equation requires that there be 2 geometric dimensions."
2074  CALL flagerror(local_error,err,error,*999)
2075  ENDIF
2076  !Set number of analytic field components
2077  number_of_analytic_components=0
2078  !Set analytic function type
2079  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE= &
2082  !Check that domain is 3D
2083  IF(number_of_dimensions/=3) THEN
2084  local_error="The number of geometric dimensions of "// &
2085  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2086  & " is invalid. The analytic function type of "// &
2087  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2088  & " for a multi-compartment diffusion equation requires that there be 3 geometric dimensions."
2089  CALL flagerror(local_error,err,error,*999)
2090  ENDIF
2091  !Set number of analytic field components
2092  number_of_analytic_components=0
2093  !Set analytic function type
2094  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE= &
2097  !Check that domain is 3D
2098  IF(number_of_dimensions/=3) THEN
2099  local_error="The number of geometric dimensions of "// &
2100  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2101  & " is invalid. The analytic function type of "// &
2102  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2103  & " for a multi-compartment diffusion equation requires that there be 3 geometric dimensions."
2104  CALL flagerror(local_error,err,error,*999)
2105  ENDIF
2106  !Set number of analytic field components
2107  number_of_analytic_components=0
2108  !Set analytic function type
2109  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE= &
2112  !Check that domain is 3D
2113  IF(number_of_dimensions/=3) THEN
2114  local_error="The number of geometric dimensions of "// &
2115  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2116  & " is invalid. The analytic function type of "// &
2117  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2118  & " for a multi-compartment diffusion requires that there be 3 geometric dimensions."
2119  CALL flagerror(local_error,err,error,*999)
2120  ENDIF
2121  !Set number of analytic field components
2122  number_of_analytic_components=0
2123  !Set analytic function type
2124  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE= &
2126  CASE DEFAULT
2127  local_error="The specified analytic function type of "// &
2128  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2129  & " is invalid for a multi-compartment diffusion equation."
2130  CALL flagerror(local_error,err,error,*999)
2131  END SELECT
2132  CASE DEFAULT
2133  local_error="The equation set subtype of "// &
2134  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2135  & " is invalid for an analytical diffusion equation."
2136  CALL flagerror(local_error,err,error,*999)
2137  END SELECT
2138  !Create analytic field if required
2139  IF(number_of_analytic_components>=1) THEN
2140  IF(equations_analytic%ANALYTIC_FIELD_AUTO_CREATED) THEN
2141  !Create the auto created source field
2142  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
2143  & equations_analytic%ANALYTIC_FIELD,err,error,*999)
2144  CALL field_label_set(equations_analytic%ANALYTIC_FIELD,"Analytic Field",err,error,*999)
2145  CALL field_type_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_general_type,err,error,*999)
2146  CALL field_dependent_type_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_independent_type, &
2147  & err,error,*999)
2148  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
2149  & err,error,*999)
2150  CALL field_mesh_decomposition_set_and_lock(equations_analytic%ANALYTIC_FIELD, &
2151  & geometric_decomposition,err,error,*999)
2152  CALL field_geometric_field_set_and_lock(equations_analytic%ANALYTIC_FIELD,equations_set%GEOMETRY% &
2153  & geometric_field,err,error,*999)
2154  CALL field_number_of_variables_set_and_lock(equations_analytic%ANALYTIC_FIELD,1,err,error,*999)
2155  CALL field_variable_types_set_and_lock(equations_analytic%ANALYTIC_FIELD,[field_u_variable_type], &
2156  & err,error,*999)
2157  CALL field_variable_label_set(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2158  & "Analytic",err,error,*999)
2159  CALL field_dimension_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2160  & field_vector_dimension_type,err,error,*999)
2161  CALL field_data_type_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2162  & field_dp_type,err,error,*999)
2163  !Set the number of analytic components
2164  CALL field_number_of_components_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2165  & number_of_analytic_components,err,error,*999)
2166  !Default the analytic components to the 1st geometric interpolation setup with constant interpolation
2167  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, &
2168  & field_u_variable_type,1,geometric_mesh_component,err,error,*999)
2169  DO component_idx=1,number_of_analytic_components
2170  CALL field_component_mesh_component_set(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2171  & component_idx,geometric_mesh_component,err,error,*999)
2172  CALL field_component_interpolation_set(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2173  & component_idx,field_constant_interpolation,err,error,*999)
2174  ENDDO !component_idx
2175  !Default the field scaling to that of the geometric field
2176  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
2177  & err,error,*999)
2178  CALL field_scaling_type_set(equations_analytic%ANALYTIC_FIELD,geometric_scaling_type,err,error,*999)
2179  ELSE
2180  !Check the user specified field
2181  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
2182  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
2183  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
2184  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
2185  IF(number_of_analytic_components==1) THEN
2186  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type, &
2187  & field_scalar_dimension_type,err,error,*999)
2188  ELSE
2189  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type, &
2190  & field_vector_dimension_type,err,error,*999)
2191  ENDIF
2192  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type, &
2193  & err,error,*999)
2194  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
2195  & number_of_analytic_components,err,error,*999)
2196  ENDIF
2197  ENDIF
2198  ELSE
2199  CALL flagerror("Equations set materials has not been finished.",err,error,*999)
2200  ENDIF
2201  ELSE
2202  CALL flagerror("Equations set materials is not associated.",err,error,*999)
2203  ENDIF
2204  ELSE
2205  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
2206  ENDIF
2207  ELSE
2208  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
2209  ENDIF
2210  ELSE
2211  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2212  ENDIF
2213  ELSE
2214  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
2215  ENDIF
2217  equations_analytic=>equations_set%ANALYTIC
2218  IF(ASSOCIATED(equations_analytic)) THEN
2219  analytic_field=>equations_analytic%ANALYTIC_FIELD
2220  IF(ASSOCIATED(analytic_field)) THEN
2221  IF(equations_analytic%ANALYTIC_FIELD_AUTO_CREATED) THEN
2222  !Finish creating the analytic field
2223  CALL field_create_finish(equations_analytic%ANALYTIC_FIELD,err,error,*999)
2224  !Set the default values for the analytic field
2225  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2226  & number_of_dimensions,err,error,*999)
2227  SELECT CASE(equations_set%SPECIFICATION(3))
2229  SELECT CASE(equations_analytic%ANALYTIC_FUNCTION_TYPE)
2231  !Set A
2232  CALL field_component_values_initialise(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2233  & field_values_set_type,1,1.0_dp,err,error,*999)
2234  !Set B
2235  CALL field_component_values_initialise(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2236  & field_values_set_type,2,1.0_dp,err,error,*999)
2237  !Set C
2238  CALL field_component_values_initialise(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2239  & field_values_set_type,3,1.0_dp,err,error,*999)
2240  !Set L
2241  CALL field_component_values_initialise(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2242  & field_values_set_type,4,1.0_dp,err,error,*999)
2244  !Do nothing
2245  CASE DEFAULT
2246  local_error="The specified analytic function type of "// &
2247  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2248  & " is invalid for a no source diffusion equation."
2249  CALL flagerror(local_error,err,error,*999)
2250  END SELECT
2252  !Do nothing
2254  !Do nothing
2256  !Do nothing
2257  CASE DEFAULT
2258  local_error="The equation set subtype of "// &
2259  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2260  & " is invalid for an analytical linear diffusion equation."
2261  CALL flagerror(local_error,err,error,*999)
2262  END SELECT
2263  ENDIF
2264  ENDIF
2265  ELSE
2266  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
2267  ENDIF
2268  CASE DEFAULT
2269  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2270  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2271  & " is invalid for a linear diffusion equation."
2272  CALL flagerror(local_error,err,error,*999)
2273  END SELECT
2274  !-----------------------------------------------------------------
2275  ! E q u a t i o n s t y p e
2276  !-----------------------------------------------------------------
2278  SELECT CASE(equations_set_setup%ACTION_TYPE)
2280  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
2281  CALL equations_create_start(equations_set,equations,err,error,*999)
2282  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
2284  ELSE
2285  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2286  ENDIF
2288  SELECT CASE(equations_set%SOLUTION_METHOD)
2290  !Finish the equations
2291  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
2292  CALL equations_create_finish(equations,err,error,*999)
2293  !Create the equations mapping.
2294  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
2295  CALL equations_mapping_dynamic_matrices_set(equations_mapping,.true.,.true.,err,error,*999)
2296  SELECT CASE(equations_set%SPECIFICATION(3))
2298  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_v_variable_type,err,error,*999)
2299  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_delvdeln_variable_type,err,error,*999)
2301  CALL flagerror("Not implemented.",err,error,*999)
2303  equations_set_field_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
2304  CALL field_parameter_set_data_get(equations_set_field_field,field_u_variable_type, &
2305  & field_values_set_type,equations_set_field_data,err,error,*999)
2306  imy_matrix = equations_set_field_data(1)
2307  ncompartments = equations_set_field_data(2)
2308  CALL equationsmapping_linearmatricesnumberset(equations_mapping,ncompartments-1,err,error,*999)
2309 
2310  ALLOCATE(variable_types(2*ncompartments))
2311  ALLOCATE(variable_u_types(ncompartments-1))
2312  DO num_var=1,ncompartments
2313  variable_types(2*num_var-1)=field_u_variable_type+(field_number_of_variable_subtypes*(num_var-1))
2314  variable_types(2*num_var)=field_deludeln_variable_type+(field_number_of_variable_subtypes*(num_var-1))
2315  ENDDO
2316  num_var_count=0
2317  DO num_var=1,ncompartments
2318  IF(num_var/=imy_matrix)THEN
2319  num_var_count=num_var_count+1
2320  variable_u_types(num_var_count)=variable_types(2*num_var-1)
2321  ENDIF
2322  ENDDO
2323  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,variable_types(2*imy_matrix-1),err,error,*999)
2324  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,variable_u_types,err,error,*999)
2325  CALL equations_mapping_rhs_variable_type_set(equations_mapping,variable_types(2*imy_matrix),err,error,*999)
2326  CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
2327  CASE DEFAULT
2328  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
2329  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
2330  END SELECT
2331  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
2332  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
2333  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
2334  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
2335  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
2336  CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
2337  ENDIF
2338  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
2339  !Create the equations matrices
2340  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
2341  !Set up matrix storage and structure
2342  IF(equations%LUMPING_TYPE==equations_lumped_matrices) THEN
2343  !Set up lumping
2344  CALL equations_matrices_dynamic_lumping_type_set(equations_matrices, &
2346  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
2348  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
2350  ELSE
2351  SELECT CASE(equations%SPARSITY_TYPE)
2353  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
2356  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
2358  & err,error,*999)
2359  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
2361  IF(equations_set%SPECIFICATION(3)==equations_set_multi_comp_transport_diffusion_subtype)THEN
2362  ALLOCATE(coupling_matrix_storage_type(ncompartments-1))
2363  ALLOCATE(coupling_matrix_structure_type(ncompartments-1))
2364  DO num_var=1,ncompartments-1
2365  coupling_matrix_storage_type(num_var)=distributed_matrix_compressed_row_storage_type
2366  coupling_matrix_structure_type(num_var)=equations_matrix_fem_structure
2367  ENDDO
2368  CALL equations_matrices_linear_storage_type_set(equations_matrices, &
2369  & coupling_matrix_storage_type, &
2370  & err,error,*999)
2371  CALL equationsmatrices_linearstructuretypeset(equations_matrices, &
2372  coupling_matrix_structure_type,err,error,*999)
2373  ENDIF
2374  CASE DEFAULT
2375  local_error="The equations matrices sparsity type of "// &
2376  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
2377  CALL flagerror(local_error,err,error,*999)
2378  END SELECT
2379  ENDIF
2380  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
2382  CALL flagerror("Not implemented.",err,error,*999)
2384  CALL flagerror("Not implemented.",err,error,*999)
2386  CALL flagerror("Not implemented.",err,error,*999)
2388  CALL flagerror("Not implemented.",err,error,*999)
2390  CALL flagerror("Not implemented.",err,error,*999)
2391  CASE DEFAULT
2392  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2393  & " is invalid."
2394  CALL flagerror(local_error,err,error,*999)
2395  END SELECT
2396  CASE DEFAULT
2397  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2398  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2399  & " is invalid for a linear diffusion equation."
2400  CALL flagerror(local_error,err,error,*999)
2401  END SELECT
2402  CASE DEFAULT
2403  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2404  & " is invalid for a linear diffusion equation."
2405  CALL flagerror(local_error,err,error,*999)
2406  END SELECT
2407  ELSE
2408  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2409  & " is not a linear diffusion equation subtype."
2410  CALL flagerror(local_error,err,error,*999)
2411  ENDIF
2412  ELSE
2413  CALL flagerror("Equations set is not associated.",err,error,*999)
2414  ENDIF
2415 
2416  exits("Diffusion_EquationsSetLinearSetup")
2417  RETURN
2418 999 errorsexits("Diffusion_EquationsSetLinearSetup",err,error)
2419  RETURN 1
2420 
2421  END SUBROUTINE diffusion_equationssetlinearsetup
2422 
2423  !
2424  !================================================================================================================================
2425  !
2426 
2428  SUBROUTINE diffusion_equationssetnonlinearsetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
2430  !Argument variables
2431  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2432  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
2433  INTEGER(INTG), INTENT(OUT) :: ERR
2434  TYPE(varying_string), INTENT(OUT) :: ERROR
2435  !Local Variables
2436  INTEGER(INTG) :: component_idx,GEOMETRIC_MESH_COMPONENT,GEOMETRIC_SCALING_TYPE,NUMBER_OF_ANALYTIC_COMPONENTS, &
2437  & NUMBER_OF_DIMENSIONS,NUMBER_OF_MATERIALS_COMPONENTS
2438  REAL(DP) :: A_PARAM,B_PARAM,C_PARAM
2439  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
2440  TYPE(equations_type), POINTER :: EQUATIONS
2441  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
2442  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
2443  TYPE(equations_set_analytic_type), POINTER :: EQUATIONS_ANALYTIC
2444  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
2445  TYPE(field_type), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,GEOMETRIC_FIELD
2446  TYPE(varying_string) :: LOCAL_ERROR
2447 
2448  enters("Diffusion_EquationsSetNonlinearSetup",err,error,*999)
2449 
2450  NULLIFY(equations)
2451  NULLIFY(equations_mapping)
2452  NULLIFY(equations_matrices)
2453  NULLIFY(geometric_decomposition)
2454 
2455  IF(ASSOCIATED(equations_set)) THEN
2456  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
2457  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
2458  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
2459  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
2460  & err,error,*999)
2461  END IF
2462  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_diffusion_subtype.OR. &
2463  & equations_set%SPECIFICATION(3)==equations_set_exponential_source_diffusion_subtype) THEN
2464  SELECT CASE(equations_set_setup%SETUP_TYPE)
2466  SELECT CASE(equations_set_setup%ACTION_TYPE)
2469  & err,error,*999)
2471  !Do nothing
2472  CASE DEFAULT
2473  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2474  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2475  & " is invalid for a nonlinear diffusion equation."
2476  CALL flagerror(local_error,err,error,*999)
2477  END SELECT
2479  !Do nothing
2480  !-----------------------------------------------------------------
2481  ! D e p e n d e n t f i e l d
2482  !-----------------------------------------------------------------
2484  SELECT CASE(equations_set_setup%ACTION_TYPE)
2486  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
2487  !Create the auto created dependent field
2488  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
2489  & dependent_field,err,error,*999)
2490  CALL field_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,"Dependent Field",err,error,*999)
2491  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
2492  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
2493  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
2494  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
2495  & err,error,*999)
2496  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
2497  & geometric_field,err,error,*999)
2498  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
2499  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,[field_u_variable_type, &
2500  & field_deludeln_variable_type],err,error,*999)
2501  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2502  & "U",err,error,*999)
2503  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
2504  & "del U/del n",err,error,*999)
2505  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2506  & field_scalar_dimension_type,err,error,*999)
2507  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
2508  & field_scalar_dimension_type,err,error,*999)
2509  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
2510  & field_dp_type,err,error,*999)
2511  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
2512  & field_dp_type,err,error,*999)
2513  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2514  & err,error,*999)
2515  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2516  & field_deludeln_variable_type,1,err,error,*999)
2517  !Default to the geometric interpolation setup
2518  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
2519  & geometric_mesh_component,err,error,*999)
2520  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
2521  & geometric_mesh_component,err,error,*999)
2522  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
2523  & geometric_mesh_component,err,error,*999)
2524  SELECT CASE(equations_set%SOLUTION_METHOD)
2526  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2527  & field_u_variable_type,1,field_node_based_interpolation,err,error,*999)
2528  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
2529  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
2530  !Default the scaling to the geometric field scaling
2531  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
2532  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
2534  CALL flagerror("Not implemented.",err,error,*999)
2536  CALL flagerror("Not implemented.",err,error,*999)
2538  CALL flagerror("Not implemented.",err,error,*999)
2540  CALL flagerror("Not implemented.",err,error,*999)
2542  CALL flagerror("Not implemented.",err,error,*999)
2543  CASE DEFAULT
2544  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2545  & " is invalid."
2546  CALL flagerror(local_error,err,error,*999)
2547  END SELECT
2548  ELSE
2549  !Check the user specified field
2550  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
2551  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
2552  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
2553  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type, &
2554  & field_deludeln_variable_type],err,error,*999)
2555  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type, &
2556  & err,error,*999)
2557  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
2558  & err,error,*999)
2559  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
2560  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
2561  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2562  & number_of_dimensions,err,error,*999)
2563  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions, &
2564  & err,error,*999)
2565  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
2566  & number_of_dimensions,err,error,*999)
2567  SELECT CASE(equations_set%SOLUTION_METHOD)
2569  DO component_idx=1,number_of_dimensions
2570  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,component_idx, &
2571  & field_node_based_interpolation,err,error,*999)
2572  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type, &
2573  & component_idx,field_node_based_interpolation,err,error,*999)
2574  ENDDO !component_idx
2576  CALL flagerror("Not implemented.",err,error,*999)
2578  CALL flagerror("Not implemented.",err,error,*999)
2580  CALL flagerror("Not implemented.",err,error,*999)
2582  CALL flagerror("Not implemented.",err,error,*999)
2584  CALL flagerror("Not implemented.",err,error,*999)
2585  CASE DEFAULT
2586  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
2587  & " is invalid."
2588  CALL flagerror(local_error,err,error,*999)
2589  END SELECT
2590  ENDIF
2592  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
2593  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
2594  ENDIF
2595  CASE DEFAULT
2596  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2597  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2598  & " is invalid for a nonlinear diffusion equation"
2599  CALL flagerror(local_error,err,error,*999)
2600  END SELECT
2601  !-----------------------------------------------------------------
2602  ! M a t e r i a l s f i e l d
2603  !-----------------------------------------------------------------
2605  SELECT CASE(equations_set_setup%ACTION_TYPE)
2607  equations_materials=>equations_set%MATERIALS
2608  IF(ASSOCIATED(equations_materials)) THEN
2609  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
2610  !Create the auto created materials field
2611  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
2612  & materials_field,err,error,*999)
2613  CALL field_label_set(equations_materials%MATERIALS_FIELD,"Materials Field",err,error,*999)
2614  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
2615  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
2616  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
2617  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
2618  & err,error,*999)
2619  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
2620  & geometric_field,err,error,*999)
2621  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
2622  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type], &
2623  & err,error,*999)
2624  CALL field_variable_label_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2625  & "Materials",err,error,*999)
2626  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2627  & field_vector_dimension_type,err,error,*999)
2628  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2629  & field_dp_type,err,error,*999)
2630  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2631  & number_of_dimensions,err,error,*999)
2632  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_diffusion_subtype) THEN
2633  !Quadratic source. Materials field components are 1 for each dimension and 3 for the quadratic source
2634  !i.e., k and a, b and c in del u/del t = div(k.grad(u(x)))+a(x)+b(x)u(x)+c(x)u^2(x)
2635  number_of_materials_components=number_of_dimensions+3
2636  ELSE
2637  !Exponential source. Matierals field components are 1 for each dimension and 3 for the exponential source
2638  !i.e., k, a, b and c in del u/del t = div(k.grad(u(x)))+a(x)+b(x)e^[c(x)u(x)]
2639  number_of_materials_components=number_of_dimensions+3
2640  ENDIF
2641  !Set the number of materials components
2642  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2643  & number_of_materials_components,err,error,*999)
2644  !Default the k materials components to the geometric interpolation setup with constant interpolation
2645  DO component_idx=1,number_of_dimensions
2646  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2647  & component_idx,geometric_mesh_component,err,error,*999)
2648  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2649  & component_idx,field_constant_interpolation,err,error,*999)
2650  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2651  & component_idx,geometric_mesh_component,err,error,*999)
2652  ENDDO !component_idx
2653  !Default the source materials components to the first component geometric interpolation with constant
2654  !interpolation
2655  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2656  & 1,geometric_mesh_component,err,error,*999)
2657  DO component_idx=number_of_dimensions+1,number_of_materials_components
2658  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2659  & component_idx,geometric_mesh_component,err,error,*999)
2660  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2661  & component_idx,field_constant_interpolation,err,error,*999)
2662  ENDDO !components_idx
2663  !Default the field scaling to that of the geometric field
2664  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
2665  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
2666  ELSE
2667  !Check the user specified field
2668  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
2669  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
2670  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
2671  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
2672  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
2673  & err,error,*999)
2674  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
2675  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2676  & number_of_dimensions,err,error,*999)
2677  number_of_materials_components=number_of_dimensions+2
2678  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
2679  & number_of_materials_components,err,error,*999)
2680  ENDIF
2681  ELSE
2682  CALL flagerror("Equations set materials is not associated.",err,error,*999)
2683  ENDIF
2685  equations_materials=>equations_set%MATERIALS
2686  IF(ASSOCIATED(equations_materials)) THEN
2687  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
2688  !Finish creating the materials field
2689  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
2690  !Set the default values for the materials field
2691  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
2692  & number_of_dimensions,err,error,*999)
2693  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_diffusion_subtype) THEN
2694  !Quadratic source. Materials field components are 1 for each dimension and 3 for the quadratic source
2695  !i.e., k and a, b and c in del u/del t = div(k.grad(u(x)))+a(x)+b(x)u(x)+c(x)u^2(x)
2696  number_of_materials_components=number_of_dimensions+3
2697  ELSE
2698  !Exponential source. Matierals field components are 1 for each dimension and 3 for the exponential source
2699  !i.e., k, a, b and c in del u/del t = div(k.grad(u(x)))+a(x)+b(x)e^[c(x)u(x)]
2700  number_of_materials_components=number_of_dimensions+3
2701  ENDIF
2702  !First set the k values to 1.0
2703  DO component_idx=1,number_of_dimensions
2704  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2705  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
2706  ENDDO !component_idx
2707  !Set the source values to 1.0
2708  DO component_idx=number_of_dimensions+1,number_of_materials_components
2709  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2710  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
2711  ENDDO !component_idx
2712  ENDIF
2713  ELSE
2714  CALL flagerror("Equations set materials is not associated.",err,error,*999)
2715  ENDIF
2716  CASE DEFAULT
2717  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2718  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2719  & " is invalid for a linear diffusion equation."
2720  CALL flagerror(local_error,err,error,*999)
2721  END SELECT
2722  !-----------------------------------------------------------------
2723  ! S o u r c e f i e l d
2724  !-----------------------------------------------------------------
2726  SELECT CASE(equations_set_setup%ACTION_TYPE)
2728  !Do nothing put the constant source directly into the RHS
2730  !Do nothing put the constant source directly into the RHS
2731  CASE DEFAULT
2732  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2733  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2734  & " is invalid for a linear diffusion equation."
2735  CALL flagerror(local_error,err,error,*999)
2736  END SELECT
2737  !-----------------------------------------------------------------
2738  ! A n a l y t i c t y p e
2739  !-----------------------------------------------------------------
2741  SELECT CASE(equations_set_setup%ACTION_TYPE)
2743  equations_analytic=>equations_set%ANALYTIC
2744  IF(ASSOCIATED(equations_analytic)) THEN
2745  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
2746  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
2747  IF(ASSOCIATED(dependent_field)) THEN
2748  equations_materials=>equations_set%MATERIALS
2749  IF(ASSOCIATED(equations_materials)) THEN
2750  IF(equations_materials%MATERIALS_FINISHED) THEN
2751  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
2752  IF(ASSOCIATED(geometric_field)) THEN
2753  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions, &
2754  & err,error,*999)
2755  IF(equations_set%SPECIFICATION(3)==equations_set_quadratic_source_diffusion_subtype) THEN
2756  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
2758  !Check that domain is 1D
2759  IF(number_of_dimensions/=1) THEN
2760  local_error="The number of geometric dimensions of "// &
2761  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2762  & " is invalid. The analytic function type of "// &
2763  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2764  & " requires that there be 1 geometric dimension."
2765  CALL flagerror(local_error,err,error,*999)
2766  ENDIF
2767  !Check the materials values are constant
2768  CALL field_component_interpolation_check(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2769  & 1,field_constant_interpolation,err,error,*999)
2770  CALL field_component_interpolation_check(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2771  & 2,field_constant_interpolation,err,error,*999)
2772  CALL field_component_interpolation_check(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2773  & 3,field_constant_interpolation,err,error,*999)
2774  !Check that the a parameter is zero.
2775  CALL field_parameter_set_get_constant(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2776  & field_values_set_type,1,a_param,err,error,*999)
2777  IF(abs(a_param)>zero_tolerance) &
2778  & CALL flagerror("The 1st material component must be zero.",err,error,*999)
2779  !Check that the b parameter is not zero.
2780  CALL field_parameter_set_get_constant(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2781  & field_values_set_type,2,b_param,err,error,*999)
2782  IF(b_param<zero_tolerance) &
2783  & CALL flagerror("The 2nd material component must be greater than zero.",err,error,*999)
2784  !Check to ensure we get real solutions
2785  CALL field_parameter_set_get_constant(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2786  & field_values_set_type,2,b_param,err,error,*999)
2787  IF((b_param*c_param)>zero_tolerance) &
2788  & CALL flagerror("The product of the 2nd and 3rd material components must not be positive.", &
2789  & err,error,*999)
2790  !Set the number of analytic field components
2791  number_of_analytic_components=1
2792  !Set analytic function type
2793  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE= &
2795  CASE DEFAULT
2796  local_error="The specified analytic function type of "// &
2797  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2798  & " is invalid for a nonlinear diffusion equation with a quadratic source."
2799  CALL flagerror(local_error,err,error,*999)
2800  END SELECT
2801  ELSE
2802  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
2804  !Check that domain is 1D
2805  IF(number_of_dimensions/=1) THEN
2806  local_error="The number of geometric dimensions of "// &
2807  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
2808  & " is invalid. The analytic function type of "// &
2809  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2810  & " requires that there be 1 geometric dimension."
2811  CALL flagerror(local_error,err,error,*999)
2812  ENDIF
2813  !Check the materials values are constant
2814  CALL field_component_interpolation_check(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2815  & 1,field_constant_interpolation,err,error,*999)
2816  CALL field_component_interpolation_check(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2817  & 2,field_constant_interpolation,err,error,*999)
2818  CALL field_component_interpolation_check(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2819  & 3,field_constant_interpolation,err,error,*999)
2820  !Check that the a parameter is not zero.
2821  CALL field_parameter_set_get_constant(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2822  & field_values_set_type,1,a_param,err,error,*999)
2823  IF(abs(a_param)<zero_tolerance) &
2824  & CALL flagerror("The 1st material component must not be zero.",err,error,*999)
2825  !Check that the c parameter is not zero.
2826  CALL field_parameter_set_get_constant(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2827  & field_values_set_type,3,c_param,err,error,*999)
2828  IF(abs(c_param)<zero_tolerance) &
2829  & CALL flagerror("The 3rd material component must not be zero.",err,error,*999)
2830  !Check to ensure we get real solutions
2831  CALL field_parameter_set_get_constant(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
2832  & field_values_set_type,2,b_param,err,error,*999)
2833  IF((a_param*b_param)>zero_tolerance) &
2834  & CALL flagerror("The product of the 1st and 2nd material components must not be positive.", &
2835  & err,error,*999)
2836  IF((a_param*c_param)<zero_tolerance) &
2837  & CALL flagerror("The product of the 1st and 3rd material components must not be negative.", &
2838  & err,error,*999)
2839  !Set the number of analytic field components
2840  number_of_analytic_components=0
2841  !Set analytic function type
2842  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE= &
2844  CASE DEFAULT
2845  local_error="The specified analytic function type of "// &
2846  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
2847  & " is invalid for a nonlinear diffusion equation with an exponential source."
2848  CALL flagerror(local_error,err,error,*999)
2849  END SELECT
2850  ENDIF
2851  !Create analytic field if required
2852  IF(number_of_analytic_components>=1) THEN
2853  IF(equations_analytic%ANALYTIC_FIELD_AUTO_CREATED) THEN
2854  !Create the auto created source field
2855  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
2856  & equations_analytic%ANALYTIC_FIELD,err,error,*999)
2857  CALL field_label_set(equations_analytic%ANALYTIC_FIELD,"Analytic Field",err,error,*999)
2858  CALL field_type_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_general_type,err,error,*999)
2859  CALL field_dependent_type_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_independent_type, &
2860  & err,error,*999)
2861  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition, &
2862  & err,error,*999)
2863  CALL field_mesh_decomposition_set_and_lock(equations_analytic%ANALYTIC_FIELD, &
2864  & geometric_decomposition,err,error,*999)
2865  CALL field_geometric_field_set_and_lock(equations_analytic%ANALYTIC_FIELD,equations_set%GEOMETRY% &
2866  & geometric_field,err,error,*999)
2867  CALL field_number_of_variables_set_and_lock(equations_analytic%ANALYTIC_FIELD,1,err,error,*999)
2868  CALL field_variable_types_set_and_lock(equations_analytic%ANALYTIC_FIELD,[field_u_variable_type], &
2869  & err,error,*999)
2870  CALL field_variable_label_set(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2871  & "Analytic",err,error,*999)
2872  CALL field_dimension_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2873  & field_vector_dimension_type,err,error,*999)
2874  CALL field_data_type_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2875  & field_dp_type,err,error,*999)
2876  !Set the number of analytic components
2877  CALL field_number_of_components_set_and_lock(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2878  & number_of_analytic_components,err,error,*999)
2879  !Default the analytic components to the 1st geometric interpolation setup with constant interpolation
2880  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD, &
2881  & field_u_variable_type,1,geometric_mesh_component,err,error,*999)
2882  DO component_idx=1,number_of_analytic_components
2883  CALL field_component_mesh_component_set(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2884  & component_idx,geometric_mesh_component,err,error,*999)
2885  CALL field_component_interpolation_set(equations_analytic%ANALYTIC_FIELD,field_u_variable_type, &
2886  & component_idx,field_constant_interpolation,err,error,*999)
2887  ENDDO !component_idx
2888  !Default the field scaling to that of the geometric field
2889  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
2890  & err,error,*999)
2891  CALL field_scaling_type_set(equations_analytic%ANALYTIC_FIELD,geometric_scaling_type,err,error,*999)
2892  ELSE
2893  !Check the user specified field
2894  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
2895  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
2896  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
2897  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
2898  IF(number_of_analytic_components==1) THEN
2899  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type, &
2900  & field_scalar_dimension_type,err,error,*999)
2901  ELSE
2902  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type, &
2903  & field_vector_dimension_type,err,error,*999)
2904  ENDIF
2905  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type, &
2906  & err,error,*999)
2907  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type, &
2908  & number_of_analytic_components,err,error,*999)
2909  ENDIF
2910  ENDIF
2911  ELSE
2912  CALL flagerror("Equations set materials is not finished.",err,error,*999)
2913  ENDIF
2914  ELSE
2915  CALL flagerror("Equations set materials is not associated.",err,error,*999)
2916  ENDIF
2917  ELSE
2918  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
2919  ENDIF
2920  ELSE
2921  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
2922  ENDIF
2923  ELSE
2924  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2925  ENDIF
2926  ELSE
2927  CALL flagerror("Equations analytic is not associated.",err,error,*999)
2928  ENDIF
2930  equations_analytic=>equations_set%ANALYTIC
2931  IF(ASSOCIATED(equations_analytic)) THEN
2932  analytic_field=>equations_analytic%ANALYTIC_FIELD
2933  IF(ASSOCIATED(analytic_field)) THEN
2934  IF(equations_analytic%ANALYTIC_FIELD_AUTO_CREATED) THEN
2935  !Finish creating the analytic field
2936  CALL field_create_finish(equations_analytic%ANALYTIC_FIELD,err,error,*999)
2937  !Set the default values for the analytic field
2938  SELECT CASE(equations_set%SPECIFICATION(3))
2940  !Do nothing
2942  !Do nothing
2943  CASE DEFAULT
2944  local_error="The equation set subtype of "// &
2945  & trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
2946  & " is invalid for an analytical nonlinear diffusion equation."
2947  CALL flagerror(local_error,err,error,*999)
2948  END SELECT
2949  ENDIF
2950  ENDIF
2951  ELSE
2952  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
2953  ENDIF
2954  CASE DEFAULT
2955  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
2956  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
2957  & " is invalid for a linear diffusion equation."
2958  CALL flagerror(local_error,err,error,*999)
2959  END SELECT
2960  !-----------------------------------------------------------------
2961  ! E q u a t i o n s t y p e
2962  !-----------------------------------------------------------------
2964  SELECT CASE(equations_set_setup%ACTION_TYPE)
2966  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
2967  CALL equations_create_start(equations_set,equations,err,error,*999)
2968  CALL equations_linearity_type_set(equations,equations_nonlinear,err,error,*999)
2970  ELSE
2971  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
2972  ENDIF
2974  SELECT CASE(equations_set%SOLUTION_METHOD)
2976  !Finish the equations
2977  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
2978  CALL equations_create_finish(equations,err,error,*999)
2979  !Create the equations mapping.
2980  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
2981  CALL equations_mapping_dynamic_matrices_set(equations_mapping,.true.,.true.,err,error,*999)
2982  CALL equationsmapping_residualvariabletypesset(equations_mapping,[field_u_variable_type],err,error,*999)
2983  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
2984  CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
2985  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
2986  !Create the equations matrices
2987  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
2988  ! Use the analytic Jacobian calculation
2990  & err,error,*999)
2991  !Set up matrix storage and structure
2992  IF(equations%LUMPING_TYPE==equations_lumped_matrices) THEN
2993  !Set up lumping
2994  CALL equations_matrices_dynamic_lumping_type_set(equations_matrices, &
2996  SELECT CASE(equations%SPARSITY_TYPE)
2998  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
3000  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
3003  & err,error,*999)
3005  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
3007  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
3009  CALL equationsmatrices_nonlinearstoragetypeset(equations_matrices, &
3012  & err,error,*999)
3013  CASE DEFAULT
3014  local_error="The equations matrices sparsity type of "// &
3015  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
3016  CALL flagerror(local_error,err,error,*999)
3017  END SELECT
3018  ELSE
3019  SELECT CASE(equations%SPARSITY_TYPE)
3021  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
3024  & err,error,*999)
3026  CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
3028  & err,error,*999)
3029  CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
3031  CALL equationsmatrices_nonlinearstoragetypeset(equations_matrices, &
3033  CALL equationsmatrices_nonlinearstructuretypeset(equations_matrices, &
3034  equations_matrix_fem_structure,err,error,*999)
3035  CASE DEFAULT
3036  local_error="The equations matrices sparsity type of "// &
3037  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
3038  CALL flagerror(local_error,err,error,*999)
3039  END SELECT
3040  ENDIF
3041  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
3043  CALL flagerror("Not implemented.",err,error,*999)
3045  CALL flagerror("Not implemented.",err,error,*999)
3047  CALL flagerror("Not implemented.",err,error,*999)
3049  CALL flagerror("Not implemented.",err,error,*999)
3051  CALL flagerror("Not implemented.",err,error,*999)
3052  CASE DEFAULT
3053  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
3054  & " is invalid."
3055  CALL flagerror(local_error,err,error,*999)
3056  END SELECT
3057  CASE DEFAULT
3058  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
3059  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
3060  & " is invalid for a nonlinear diffusion equation."
3061  CALL flagerror(local_error,err,error,*999)
3062  END SELECT
3063  CASE DEFAULT
3064  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
3065  & " is invalid for a nonlinear diffusion equation."
3066  CALL flagerror(local_error,err,error,*999)
3067  END SELECT
3068  ELSE
3069  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
3070  & " is not a nonlinear diffusion equation subtype."
3071  CALL flagerror(local_error,err,error,*999)
3072  ENDIF
3073  ELSE
3074  CALL flagerror("Equations set is not associated.",err,error,*999)
3075  ENDIF
3076 
3077  exits("Diffusion_EquationsSetNonlinearSetup")
3078  RETURN
3079 999 errors("Diffusion_EquationsSetNonlinearSetup",err,error)
3080  exits("Diffusion_EquationsSetNonlinearSetup")
3081  RETURN 1
3082 
3084 
3085  !
3086  !================================================================================================================================
3087  !
3088 
3090  SUBROUTINE diffusion_equation_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
3092  !Argument variables
3093  TYPE(problem_type), POINTER :: PROBLEM
3094  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
3095  INTEGER(INTG), INTENT(OUT) :: ERR
3096  TYPE(varying_string), INTENT(OUT) :: ERROR
3097  !Local Variables
3098  TYPE(varying_string) :: LOCAL_ERROR
3099 
3100  enters("DIFFUSION_EQUATION_PROBLEM_SETUP",err,error,*999)
3101 
3102  IF(ASSOCIATED(problem)) THEN
3103  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
3104  CALL flagerror("Problem specification is not allocated.",err,error,*999)
3105  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
3106  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
3107  END IF
3108  SELECT CASE(problem%SPECIFICATION(3))
3110  CALL diffusion_equation_problem_linear_setup(problem,problem_setup,err,error,*999)
3112  CALL diffusion_equation_problem_linear_setup(problem,problem_setup,err,error,*999)
3114  CALL diffusion_equation_problem_nonlinear_setup(problem,problem_setup,err,error,*999)
3116  CALL diffusion_equation_problem_linear_setup(problem,problem_setup,err,error,*999)
3118  CALL diffusion_equation_problem_linear_setup(problem,problem_setup,err,error,*999)
3120  CALL flagerror("Not implemented.",err,error,*999)
3121  CASE DEFAULT
3122  local_error="Problem subtype "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
3123  & " is not valid for a diffusion equation type of a classical field problem class."
3124  CALL flagerror(local_error,err,error,*999)
3125  END SELECT
3126  ELSE
3127  CALL flagerror("Problem is not associated.",err,error,*999)
3128  ENDIF
3129 
3130  exits("DIFFUSION_EQUATION_PROBLEM_SETUP")
3131  RETURN
3132 999 errorsexits("DIFFUSION_EQUATION_PROBLEM_SETUP",err,error)
3133  RETURN 1
3134  END SUBROUTINE diffusion_equation_problem_setup
3135 
3136  !
3137  !================================================================================================================================
3138  !
3139 
3141  SUBROUTINE diffusion_equation_pre_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
3143  !Argument variables
3144  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
3145  TYPE(solver_type), POINTER :: SOLVER
3146  INTEGER(INTG), INTENT(OUT) :: ERR
3147  TYPE(varying_string), INTENT(OUT) :: ERROR
3148 
3149  !Local variables
3150  LOGICAL :: UPDATE_MATERIALS
3151  LOGICAL :: UPDATE_BOUNDARY_CONDITIONS
3152  TYPE(varying_string) :: LOCAL_ERROR
3153 
3154  update_materials = .false.
3155  update_boundary_conditions = .true.
3156 
3157  IF( update_materials ) THEN
3158  !CALL DIFFUSION_EQUATION_PRE_SOLVE_UPDATE_MATERIALS_FIELD(CONTROL_LOOP,SOLVER,ERR,ERROR,*999)
3159  ENDIF
3160 
3161  ! IF( UPDATE_BOUNDARY_CONDITIONS ) THEN
3162  ! CALL Diffusion_PreSolveUpdateBoundaryConditions(CONTROL_LOOP,SOLVER,ERR,ERROR,*999)
3163  ! ENDIF
3164 
3165  IF(ASSOCIATED(control_loop)) THEN
3166  IF(ASSOCIATED(solver)) THEN
3167  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
3168  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
3169  CALL flagerror("Problem specification is not allocated.",err,error,*999)
3170  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
3171  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
3172  END IF
3173  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
3176  ! do nothing ???
3177  CALL diffusion_presolveupdateanalyticvalues(control_loop,solver,err,error,*999)
3180  CALL write_string(general_output_type,"ALE diffusion pre solve... ",err,error,*999)
3181  IF(solver%DYNAMIC_SOLVER%ALE) THEN
3182  !First update mesh and calculate boundary velocity values
3183  CALL diffusion_presolvealeupdatemesh(control_loop,solver,err,error,*999)
3184  !Then apply both normal and moving mesh boundary conditions
3185  !CALL DIFFUSION_PRE_SOLVE_UPDATE_BOUNDARY_CONDITIONS(CONTROL_LOOP,SOLVER,ERR,ERROR,*999)
3186  ELSE
3187  CALL flagerror("Mesh motion calculation not successful for ALE problem.",err,error,*999)
3188  END IF
3189  CASE DEFAULT
3190  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
3191  & " is not valid for a diffusion equation type of a classical field problem class."
3192  CALL flagerror(local_error,err,error,*999)
3193  END SELECT
3194  ELSE
3195  CALL flagerror("Problem is not associated.",err,error,*999)
3196  ENDIF
3197  ELSE
3198  CALL flagerror("Solver is not associated.",err,error,*999)
3199  ENDIF
3200  ELSE
3201  CALL flagerror("Control loop is not associated.",err,error,*999)
3202  ENDIF
3203 
3204  exits("DIFFUSION_EQUATION_PRE_SOLVE")
3205  RETURN
3206 999 errorsexits("DIFFUSION_EQUATION_PRE_SOLVE",err,error)
3207  RETURN 1
3208  END SUBROUTINE diffusion_equation_pre_solve
3209 
3210  !
3211  !================================================================================================================================
3212  !
3214  SUBROUTINE diffusion_presolveupdateboundaryconditions(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
3216  !Argument variables
3217  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
3218  TYPE(solver_type), POINTER :: SOLVER
3219  INTEGER(INTG), INTENT(OUT) :: ERR
3220  TYPE(varying_string), INTENT(OUT) :: ERROR
3221  !Local Variables
3222 ! TYPE(FIELD_TYPE), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
3223 ! ! TYPE(FIELD_TYPE), POINTER :: FIELD !<A pointer to the field
3224 ! TYPE(FIELD_VARIABLE_TYPE), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
3225 ! TYPE(SOLVER_EQUATIONS_TYPE), POINTER :: SOLVER_EQUATIONS !<A pointer to the solver equations
3226 ! TYPE(SOLVER_MAPPING_TYPE), POINTER :: SOLVER_MAPPING !<A pointer to the solver mapping
3227 ! TYPE(EQUATIONS_SET_TYPE), POINTER :: EQUATIONS_SET !<A pointer to the equations set
3228 ! TYPE(EQUATIONS_TYPE), POINTER :: EQUATIONS
3229 ! TYPE(DOMAIN_TYPE), POINTER :: DOMAIN
3230 ! TYPE(DOMAIN_NODES_TYPE), POINTER :: DOMAIN_NODES
3231 ! ! TYPE(DOMAIN_TOPOLOGY_TYPE), POINTER :: DOMAIN_TOPOLOGY
3232 ! TYPE(VARYING_STRING) :: LOCAL_ERROR
3233 ! ! TYPE(BOUNDARY_CONDITIONS_VARIABLE_TYPE), POINTER :: BOUNDARY_CONDITIONS_VARIABLE
3234 ! ! TYPE(BOUNDARY_CONDITIONS_TYPE), POINTER :: BOUNDARY_CONDITIONS
3235 ! ! REAL(DP), POINTER :: BOUNDARY_VALUES(:)
3236 ! REAL(DP), POINTER :: GEOMETRIC_PARAMETERS(:)
3237 ! INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,BOUNDARY_CONDITION_CHECK_VARIABLE
3238 !
3239 ! REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
3240 ! REAL(DP) :: VALUE,X(3) !<The value to add
3241 ! ! REAL(DP) :: k_xx, k_yy, k_zz
3242 ! INTEGER(INTG) :: component_idx,deriv_idx,dim_idx,local_ny,node_idx,variable_idx
3243 ! INTEGER(INTG) :: VARIABLE_TYPE !<The field variable type to add \see FIELD_ROUTINES_VariableTypes,FIELD_ROUTINES
3244 ! INTEGER(INTG) :: ANALYTIC_FUNCTION_TYPE
3245 ! INTEGER(INTG) :: GLOBAL_DERIV_INDEX
3246 ! ! INTEGER(INTG) :: FIELD_SET_TYPE !<The field parameter set identifier \see FIELD_ROUTINES_ParameterSetTypes,FIELD_ROUTINES
3247 ! ! INTEGER(INTG) :: DERIVATIVE_NUMBER !<The node derivative number
3248 ! ! INTEGER(INTG) :: COMPONENT_NUMBER !<The field variable component number
3249 ! ! INTEGER(INTG) :: TOTAL_NUMBER_OF_NODES !<The total number of (geometry) nodes
3250 ! ! INTEGER(INTG) :: LOCAL_NODE_NUMBER
3251 ! ! INTEGER(INTG) :: EQUATIONS_SET_IDX
3252 ! ! INTEGER(INTG) :: equations_row_number
3253 
3254  enters("Diffusion_PreSolveUpdateBoundaryConditions",err,error,*999)
3255 
3256  CALL flagerror("Not implemented.",err,error,*999)!This routine previously set analytic BCs, but this has been moved. Needs rewriting to set
3257  !boundary conditions from file, time varying if appropriate.
3258 
3259 ! IF(ASSOCIATED(CONTROL_LOOP)) THEN
3260 ! CALL CONTROL_LOOP_CURRENT_TIMES_GET(CONTROL_LOOP,CURRENT_TIME,TIME_INCREMENT,ERR,ERROR,*999)
3261 ! !write(*,*)'CURRENT_TIME = ',CURRENT_TIME
3262 ! !write(*,*)'TIME_INCREMENT = ',TIME_INCREMENT
3263 ! IF(ASSOCIATED(SOLVER)) THEN
3264 ! IF(ASSOCIATED(CONTROL_LOOP%PROBLEM)) THEN
3265 ! SELECT CASE(CONTROL_LOOP%PROBLEM%SPECIFICATION(3))
3266 ! CASE(PROBLEM_LINEAR_SOURCE_DIFFUSION_SUBTYPE)
3267 ! SOLVER_EQUATIONS=>SOLVER%SOLVER_EQUATIONS
3268 ! IF(ASSOCIATED(SOLVER_EQUATIONS)) THEN
3269 ! SOLVER_MAPPING=>SOLVER_EQUATIONS%SOLVER_MAPPING
3270 ! EQUATIONS=>SOLVER_MAPPING%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
3271 ! IF(ASSOCIATED(EQUATIONS)) THEN
3272 ! EQUATIONS_SET=>EQUATIONS%EQUATIONS_SET
3273 ! IF(ASSOCIATED(EQUATIONS_SET)) THEN
3274 ! IF(ASSOCIATED(EQUATIONS_SET%ANALYTIC)) THEN
3275 ! DEPENDENT_FIELD=>EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD
3276 ! IF(ASSOCIATED(DEPENDENT_FIELD)) THEN
3277 ! GEOMETRIC_FIELD=>EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD
3278 ! IF(ASSOCIATED(GEOMETRIC_FIELD)) THEN
3279 ! CALL FIELD_NUMBER_OF_COMPONENTS_GET(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,&
3280 ! & NUMBER_OF_DIMENSIONS,ERR,ERROR,*999)
3281 ! NULLIFY(GEOMETRIC_VARIABLE)
3282 ! NULLIFY(GEOMETRIC_PARAMETERS)
3283 ! CALL FIELD_VARIABLE_GET(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,GEOMETRIC_VARIABLE,ERR,ERROR,*999)
3284 ! CALL FIELD_PARAMETER_SET_DATA_GET(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE,&
3285 ! & GEOMETRIC_PARAMETERS,ERR,ERROR,*999)
3286 ! DO variable_idx=1,DEPENDENT_FIELD%NUMBER_OF_VARIABLES
3287 ! variable_type=DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
3288 ! FIELD_VARIABLE=>DEPENDENT_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
3289 ! IF(ASSOCIATED(FIELD_VARIABLE)) THEN
3290 ! DO component_idx=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
3291 ! IF(FIELD_VARIABLE%COMPONENTS(component_idx)%INTERPOLATION_TYPE== &
3292 ! & FIELD_NODE_BASED_INTERPOLATION) THEN
3293 ! DOMAIN=>FIELD_VARIABLE%COMPONENTS(component_idx)%DOMAIN
3294 ! IF(ASSOCIATED(DOMAIN)) THEN
3295 ! IF(ASSOCIATED(DOMAIN%TOPOLOGY)) THEN
3296 ! DOMAIN_NODES=>DOMAIN%TOPOLOGY%NODES
3297 ! IF(ASSOCIATED(DOMAIN_NODES)) THEN
3298 ! !Loop over the local nodes excluding the ghosts.
3299 ! DO node_idx=1,DOMAIN_NODES%NUMBER_OF_NODES
3300 ! !!TODO \todo We should interpolate the geometric field here and the node position.
3301 ! DO dim_idx=1,NUMBER_OF_DIMENSIONS
3302 ! local_ny= &
3303 ! & GEOMETRIC_VARIABLE%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP(1,node_idx)
3304 ! X(dim_idx)=GEOMETRIC_PARAMETERS(local_ny)
3305 ! ENDDO !dim_idx
3306 ! !Loop over the derivatives
3307 ! DO deriv_idx=1,DOMAIN_NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES
3308 ! ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE
3309 ! GLOBAL_DERIV_INDEX=DOMAIN_NODES%NODES(node_idx)%GLOBAL_DERIVATIVE_INDEX(deriv_idx)
3310 ! CALL DIFFUSION_EQUATION_ANALYTIC_FUNCTIONS(VALUE,X, &
3311 ! & CURRENT_TIME,variable_type,GLOBAL_DERIV_INDEX, &
3312 ! & ANALYTIC_FUNCTION_TYPE,ERR,ERROR,*999)
3313 ! local_ny=FIELD_VARIABLE%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
3314 ! & NODE_PARAM2DOF_MAP(deriv_idx,node_idx)
3315 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_DOF(DEPENDENT_FIELD,variable_type, &
3316 ! & FIELD_ANALYTIC_VALUES_SET_TYPE,local_ny,VALUE,ERR,ERROR,*999)
3317 ! BOUNDARY_CONDITION_CHECK_VARIABLE=SOLVER_EQUATIONS%BOUNDARY_CONDITIONS% &
3318 ! & BOUNDARY_CONDITIONS_VARIABLE_TYPE_MAP(FIELD_U_VARIABLE_TYPE)%PTR% &
3319 ! & CONDITION_TYPES(local_ny)
3320 ! IF(BOUNDARY_CONDITION_CHECK_VARIABLE==BOUNDARY_CONDITION_FIXED) THEN
3321 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_DOF(DEPENDENT_FIELD, &
3322 ! & variable_type,FIELD_VALUES_SET_TYPE,local_ny, &
3323 ! & VALUE,ERR,ERROR,*999)
3324 ! ENDIF
3325 ! ! IF(variable_type==FIELD_U_VARIABLE_TYPE) THEN
3326 ! ! IF(DOMAIN_NODES%NODES(node_idx)%BOUNDARY_NODE) THEN
3327 ! !If we are a boundary node then set the analytic value on the boundary
3328 ! ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3329 ! ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3330 ! ! ENDIF
3331 ! ! ENDIF
3332 ! ENDDO !deriv_idx
3333 ! ENDDO !node_idx
3334 ! ELSE
3335 ! CALL FlagError("Domain topology nodes is not associated.",ERR,ERROR,*999)
3336 ! ENDIF
3337 ! ELSE
3338 ! CALL FlagError("Domain topology is not associated.",ERR,ERROR,*999)
3339 ! ENDIF
3340 ! ELSE
3341 ! CALL FlagError("Domain is not associated.",ERR,ERROR,*999)
3342 ! ENDIF
3343 ! ELSE
3344 ! CALL FlagError("Only node based interpolation is implemented.",ERR,ERROR,*999)
3345 ! ENDIF
3346 ! ENDDO !component_idx
3347 ! CALL FIELD_PARAMETER_SET_UPDATE_START(DEPENDENT_FIELD,variable_type, &
3348 ! & FIELD_ANALYTIC_VALUES_SET_TYPE,ERR,ERROR,*999)
3349 ! CALL FIELD_PARAMETER_SET_UPDATE_FINISH(DEPENDENT_FIELD,variable_type, &
3350 ! & FIELD_ANALYTIC_VALUES_SET_TYPE,ERR,ERROR,*999)
3351 ! CALL FIELD_PARAMETER_SET_UPDATE_START(DEPENDENT_FIELD,variable_type, &
3352 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3353 ! CALL FIELD_PARAMETER_SET_UPDATE_FINISH(DEPENDENT_FIELD,variable_type, &
3354 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3355 ! ELSE
3356 ! CALL FlagError("Field variable is not associated.",ERR,ERROR,*999)
3357 ! ENDIF
3358 !
3359 ! ENDDO !variable_idx
3360 ! CALL FIELD_PARAMETER_SET_DATA_RESTORE(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,&
3361 ! & FIELD_VALUES_SET_TYPE,GEOMETRIC_PARAMETERS,ERR,ERROR,*999)
3362 ! ELSE
3363 ! CALL FlagError("Equations set geometric field is not associated.",ERR,ERROR,*999)
3364 ! ENDIF
3365 ! ELSE
3366 ! CALL FlagError("Equations set dependent field is not associated.",ERR,ERROR,*999)
3367 ! ENDIF
3368 ! ELSE
3369 ! !CALL FlagError("Equations set analytic is not associated.",ERR,ERROR,*999)
3370 ! ENDIF
3371 ! ELSE
3372 ! CALL FlagError("Equations set is not associated.",ERR,ERROR,*999)
3373 ! ENDIF
3374 ! ELSE
3375 ! CALL FlagError("Equations are not associated.",ERR,ERROR,*999)
3376 ! END IF
3377 ! ELSE
3378 ! CALL FlagError("Solver equations are not associated.",ERR,ERROR,*999)
3379 ! END IF
3380 ! CALL FIELD_PARAMETER_SET_UPDATE_START(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, &
3381 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3382 ! CALL FIELD_PARAMETER_SET_UPDATE_FINISH(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, &
3383 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3384 ! !do nothing?!
3385 ! CASE(PROBLEM_NONLINEAR_SOURCE_DIFFUSION_SUBTYPE)
3386 ! !do nothing?!
3387 ! CASE(PROBLEM_NO_SOURCE_DIFFUSION_SUBTYPE)
3388 ! SOLVER_EQUATIONS=>SOLVER%SOLVER_EQUATIONS
3389 ! IF(ASSOCIATED(SOLVER_EQUATIONS)) THEN
3390 ! SOLVER_MAPPING=>SOLVER_EQUATIONS%SOLVER_MAPPING
3391 ! EQUATIONS=>SOLVER_MAPPING%EQUATIONS_SET_TO_SOLVER_MAP(1)%EQUATIONS
3392 ! IF(ASSOCIATED(EQUATIONS)) THEN
3393 ! EQUATIONS_SET=>EQUATIONS%EQUATIONS_SET
3394 ! IF(ASSOCIATED(EQUATIONS_SET)) THEN
3395 ! IF(ASSOCIATED(EQUATIONS_SET%ANALYTIC)) THEN
3396 ! DEPENDENT_FIELD=>EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD
3397 ! IF(ASSOCIATED(DEPENDENT_FIELD)) THEN
3398 ! GEOMETRIC_FIELD=>EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD
3399 ! IF(ASSOCIATED(GEOMETRIC_FIELD)) THEN
3400 ! CALL FIELD_NUMBER_OF_COMPONENTS_GET(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,&
3401 ! & NUMBER_OF_DIMENSIONS,ERR,ERROR,*999)
3402 ! NULLIFY(GEOMETRIC_VARIABLE)
3403 ! NULLIFY(GEOMETRIC_PARAMETERS)
3404 ! CALL FIELD_VARIABLE_GET(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,GEOMETRIC_VARIABLE,ERR,ERROR,*999)
3405 ! CALL FIELD_PARAMETER_SET_DATA_GET(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE,&
3406 ! & GEOMETRIC_PARAMETERS,ERR,ERROR,*999)
3407 ! DO variable_idx=1,DEPENDENT_FIELD%NUMBER_OF_VARIABLES
3408 ! variable_type=DEPENDENT_FIELD%VARIABLES(variable_idx)%VARIABLE_TYPE
3409 ! FIELD_VARIABLE=>DEPENDENT_FIELD%VARIABLE_TYPE_MAP(variable_type)%PTR
3410 ! IF(ASSOCIATED(FIELD_VARIABLE)) THEN
3411 ! DO component_idx=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
3412 ! IF(FIELD_VARIABLE%COMPONENTS(component_idx)%INTERPOLATION_TYPE== &
3413 ! & FIELD_NODE_BASED_INTERPOLATION) THEN
3414 ! DOMAIN=>FIELD_VARIABLE%COMPONENTS(component_idx)%DOMAIN
3415 ! IF(ASSOCIATED(DOMAIN)) THEN
3416 ! IF(ASSOCIATED(DOMAIN%TOPOLOGY)) THEN
3417 ! DOMAIN_NODES=>DOMAIN%TOPOLOGY%NODES
3418 ! IF(ASSOCIATED(DOMAIN_NODES)) THEN
3419 ! !Loop over the local nodes excluding the ghosts.
3420 ! DO node_idx=1,DOMAIN_NODES%NUMBER_OF_NODES
3421 ! !!TODO \todo We should interpolate the geometric field here and the node position.
3422 ! DO dim_idx=1,NUMBER_OF_DIMENSIONS
3423 ! local_ny= &
3424 ! & GEOMETRIC_VARIABLE%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP(1,node_idx)
3425 ! X(dim_idx)=GEOMETRIC_PARAMETERS(local_ny)
3426 ! ENDDO !dim_idx
3427 ! !Loop over the derivatives
3428 ! DO deriv_idx=1,DOMAIN_NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES
3429 ! ANALYTIC_FUNCTION_TYPE=EQUATIONS_SET%ANALYTIC%ANALYTIC_FUNCTION_TYPE
3430 ! GLOBAL_DERIV_INDEX=DOMAIN_NODES%NODES(node_idx)%GLOBAL_DERIVATIVE_INDEX(deriv_idx)
3431 ! CALL DIFFUSION_EQUATION_ANALYTIC_FUNCTIONS(VALUE,X, &
3432 ! & CURRENT_TIME,variable_type,GLOBAL_DERIV_INDEX, &
3433 ! & ANALYTIC_FUNCTION_TYPE,ERR,ERROR,*999)
3434 ! local_ny=FIELD_VARIABLE%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
3435 ! & NODE_PARAM2DOF_MAP(deriv_idx,node_idx)
3436 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_DOF(DEPENDENT_FIELD,variable_type, &
3437 ! & FIELD_ANALYTIC_VALUES_SET_TYPE,local_ny,VALUE,ERR,ERROR,*999)
3438 ! BOUNDARY_CONDITION_CHECK_VARIABLE=SOLVER_EQUATIONS%BOUNDARY_CONDITIONS% &
3439 ! & BOUNDARY_CONDITIONS_VARIABLE_TYPE_MAP(FIELD_U_VARIABLE_TYPE)%PTR% &
3440 ! & CONDITION_TYPES(local_ny)
3441 ! IF(BOUNDARY_CONDITION_CHECK_VARIABLE==BOUNDARY_CONDITION_FIXED) THEN
3442 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_DOF(DEPENDENT_FIELD, &
3443 ! & variable_type,FIELD_VALUES_SET_TYPE,local_ny, &
3444 ! & VALUE,ERR,ERROR,*999)
3445 ! ENDIF
3446 ! ! IF(variable_type==FIELD_U_VARIABLE_TYPE .OR. variable_type==FIELD_V_VARIABLE_TYPE) THEN
3447 ! ! IF(DOMAIN_NODES%NODES(node_idx)%BOUNDARY_NODE) THEN
3448 ! ! If we are a boundary node then set the analytic value on the boundary
3449 ! ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3450 ! ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3451 ! ! ENDIF
3452 ! ! ENDIF
3453 ! ENDDO !deriv_idx
3454 ! ENDDO !node_idx
3455 ! ELSE
3456 ! CALL FlagError("Domain topology nodes is not associated.",ERR,ERROR,*999)
3457 ! ENDIF
3458 ! ELSE
3459 ! CALL FlagError("Domain topology is not associated.",ERR,ERROR,*999)
3460 ! ENDIF
3461 ! ELSE
3462 ! CALL FlagError("Domain is not associated.",ERR,ERROR,*999)
3463 ! ENDIF
3464 ! ELSE
3465 ! CALL FlagError("Only node based interpolation is implemented.",ERR,ERROR,*999)
3466 ! ENDIF
3467 ! ENDDO !component_idx
3468 ! CALL FIELD_PARAMETER_SET_UPDATE_START(DEPENDENT_FIELD,variable_type, &
3469 ! & FIELD_ANALYTIC_VALUES_SET_TYPE,ERR,ERROR,*999)
3470 ! CALL FIELD_PARAMETER_SET_UPDATE_FINISH(DEPENDENT_FIELD,variable_type, &
3471 ! & FIELD_ANALYTIC_VALUES_SET_TYPE,ERR,ERROR,*999)
3472 ! CALL FIELD_PARAMETER_SET_UPDATE_START(DEPENDENT_FIELD,variable_type, &
3473 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3474 ! CALL FIELD_PARAMETER_SET_UPDATE_FINISH(DEPENDENT_FIELD,variable_type, &
3475 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3476 ! ELSE
3477 ! CALL FlagError("Field variable is not associated.",ERR,ERROR,*999)
3478 ! ENDIF
3479 !
3480 ! ENDDO !variable_idx
3481 ! CALL FIELD_PARAMETER_SET_DATA_RESTORE(GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,&
3482 ! & FIELD_VALUES_SET_TYPE,GEOMETRIC_PARAMETERS,ERR,ERROR,*999)
3483 ! ELSE
3484 ! CALL FlagError("Equations set geometric field is not associated.",ERR,ERROR,*999)
3485 ! ENDIF
3486 ! ELSE
3487 ! CALL FlagError("Equations set dependent field is not associated.",ERR,ERROR,*999)
3488 ! ENDIF
3489 ! ELSE
3490 ! !CALL FlagError("Equations set analytic is not associated.",ERR,ERROR,*999)
3491 ! ENDIF
3492 ! ELSE
3493 ! CALL FlagError("Equations set is not associated.",ERR,ERROR,*999)
3494 ! ENDIF
3495 ! ELSE
3496 ! CALL FlagError("Equations are not associated.",ERR,ERROR,*999)
3497 ! END IF
3498 ! ELSE
3499 ! CALL FlagError("Solver equations are not associated.",ERR,ERROR,*999)
3500 ! END IF
3501 ! CALL FIELD_PARAMETER_SET_UPDATE_START(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, &
3502 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3503 ! CALL FIELD_PARAMETER_SET_UPDATE_FINISH(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, &
3504 ! & FIELD_VALUES_SET_TYPE,ERR,ERROR,*999)
3505 ! CASE(PROBLEM_NO_SOURCE_ALE_DIFFUSION_SUBTYPE)
3506 ! ! do nothing ???
3507 ! CASE(PROBLEM_LINEAR_SOURCE_ALE_DIFFUSION_SUBTYPE)
3508 ! ! do nothing ???
3509 ! CASE(PROBLEM_NONLINEAR_SOURCE_ALE_DIFFUSION_SUBTYPE)
3510 ! ! do nothing ???
3511 ! CASE DEFAULT
3512 ! LOCAL_ERROR="Problem subtype "//TRIM(NUMBER_TO_VSTRING(CONTROL_LOOP%PROBLEM%SPECIFICATION(3),"*",ERR,ERROR))// &
3513 ! & " is not valid for a diffusion equation type of a classical field problem class."
3514 ! CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
3515 ! END SELECT
3516 ! ELSE
3517 ! CALL FlagError("Problem is not associated.",ERR,ERROR,*999)
3518 ! ENDIF
3519 ! ELSE
3520 ! CALL FlagError("Solver is not associated.",ERR,ERROR,*999)
3521 ! ENDIF
3522 ! ELSE
3523 ! CALL FlagError("Control loop is not associated.",ERR,ERROR,*999)
3524  ! ENDIF
3525 
3526  exits("Diffusion_PreSolveUpdateBoundaryConditions")
3527  RETURN
3528 999 errors("Diffusion_PreSolveUpdateBoundaryConditions",err,error)
3529  exits("Diffusion_PreSolveUpdateBoundaryConditions")
3530  RETURN 1
3531 
3533 
3534  !
3535  !================================================================================================================================
3536  !
3537  !updates the boundary conditions and source term to the required analytic values
3538  SUBROUTINE diffusion_presolveupdateanalyticvalues(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
3540  !Argument variables
3541  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
3542  TYPE(solver_type), POINTER :: SOLVER
3543  INTEGER(INTG), INTENT(OUT) :: ERR
3544  TYPE(varying_string), INTENT(OUT) :: ERROR
3545  !Local Variables
3546  TYPE(field_type), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD,SOURCE_FIELD
3547 ! TYPE(FIELD_TYPE), POINTER :: FIELD !<A pointer to the field
3548  TYPE(field_variable_type), POINTER :: ANALYTIC_VARIABLE,FIELD_VARIABLE,GEOMETRIC_VARIABLE,MATERIALS_VARIABLE
3549  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
3550  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
3551  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
3552  TYPE(equations_set_source_type), POINTER :: EQUATIONS_SOURCE
3553  TYPE(equations_type), POINTER :: EQUATIONS
3554  TYPE(domain_type), POINTER :: DOMAIN
3555  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
3556  ! TYPE(DOMAIN_TOPOLOGY_TYPE), POINTER :: DOMAIN_TOPOLOGY
3557  TYPE(varying_string) :: LOCAL_ERROR
3558  TYPE(boundary_conditions_variable_type), POINTER :: BOUNDARY_CONDITIONS_VARIABLE
3559  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
3560 ! REAL(DP), POINTER :: BOUNDARY_VALUES(:)
3561  REAL(DP), POINTER :: ANALYTIC_PARAMETERS(:),GEOMETRIC_PARAMETERS(:),MATERIALS_PARAMETERS(:)
3562  INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,BOUNDARY_CONDITION_CHECK_VARIABLE
3563 
3564  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
3565  REAL(DP) :: NORMAL(3),TANGENTS(3,3),VALUE,X(3),VALUE_SOURCE
3566 ! REAL(DP) :: k_xx, k_yy, k_zz
3567  INTEGER(INTG) :: component_idx,deriv_idx,dim_idx,local_ny,node_idx,eqnset_idx
3568  INTEGER(INTG) :: VARIABLE_TYPE
3569  INTEGER(INTG) :: ANALYTIC_FUNCTION_TYPE
3570  INTEGER(INTG) :: GLOBAL_DERIV_INDEX
3571  REAL(DP) :: A1,D1
3572 ! INTEGER(INTG) :: FIELD_SET_TYPE !<The field parameter set identifier \see FIELD_ROUTINES_ParameterSetTypes,FIELD_ROUTINES
3573 ! INTEGER(INTG) :: DERIVATIVE_NUMBER !<The node derivative number
3574 ! INTEGER(INTG) :: COMPONENT_NUMBER !<The field variable component number
3575 ! INTEGER(INTG) :: TOTAL_NUMBER_OF_NODES !<The total number of (geometry) nodes
3576 ! INTEGER(INTG) :: LOCAL_NODE_NUMBER
3577 ! INTEGER(INTG) :: EQUATIONS_SET_IDX
3578 ! INTEGER(INTG) :: equations_row_number
3579 
3580  enters("Diffusion_PreSolveUpdateAnalyticValues",err,error,*999)
3581 
3582  a1 = 0.4_dp
3583  d1=1.0_dp
3584 
3585  IF(ASSOCIATED(control_loop)) THEN
3586  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
3587  IF(ASSOCIATED(solver)) THEN
3588  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
3589  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
3590  CALL flagerror("Problem specification is not allocated.",err,error,*999)
3591  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
3592  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
3593  END IF
3594  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
3597  solver_equations=>solver%SOLVER_EQUATIONS
3598  IF(ASSOCIATED(solver_equations)) THEN
3599  !loop over all the equation sets and set the appropriate field variable type BCs and
3600  !the source field associated with each equation set
3601  DO eqnset_idx=1,solver_equations%SOLVER_MAPPING%NUMBER_OF_EQUATIONS_SETS
3602  solver_mapping=>solver_equations%SOLVER_MAPPING
3603  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(eqnset_idx)%EQUATIONS
3604  IF(ASSOCIATED(equations)) THEN
3605  equations_set=>equations%EQUATIONS_SET
3606  IF(ASSOCIATED(equations_set)) THEN
3607  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
3608  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
3609  IF(ASSOCIATED(dependent_field)) THEN
3610  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
3611  IF(ASSOCIATED(geometric_field)) THEN
3612  analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
3613  CALL field_number_of_components_get(geometric_field,field_u_variable_type,&
3614  & number_of_dimensions,err,error,*999)
3615  NULLIFY(geometric_variable)
3616  NULLIFY(geometric_parameters)
3617  CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
3618  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,&
3619  & geometric_parameters,err,error,*999)
3620  equations_set%ANALYTIC%ANALYTIC_USER_PARAMS(1)=current_time
3621  NULLIFY(analytic_variable)
3622  NULLIFY(analytic_parameters)
3623  IF(ASSOCIATED(analytic_field)) THEN
3624  CALL field_variable_get(analytic_field,field_u_variable_type,analytic_variable,err,error,*999)
3625  CALL field_parameter_set_data_get(analytic_field,field_u_variable_type,field_values_set_type, &
3626  & analytic_parameters,err,error,*999)
3627  ENDIF
3628  NULLIFY(materials_field)
3629  NULLIFY(materials_variable)
3630  NULLIFY(materials_parameters)
3631  IF(ASSOCIATED(equations_set%MATERIALS)) THEN
3632  materials_field=>equations_set%MATERIALS%MATERIALS_FIELD
3633  CALL field_variable_get(materials_field,field_u_variable_type,materials_variable,err,error,*999)
3634  CALL field_parameter_set_data_get(materials_field,field_u_variable_type,field_values_set_type, &
3635  & materials_parameters,err,error,*999)
3636  ENDIF
3637  ! DO variable_idx=1,DEPENDENT_FIELD%NUMBER_OF_VARIABLES
3638  variable_type=dependent_field%VARIABLES(2*eqnset_idx-1)%VARIABLE_TYPE
3639  field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
3640  IF(ASSOCIATED(field_variable)) THEN
3641  boundary_conditions=>solver_equations%BOUNDARY_CONDITIONS
3642  IF(ASSOCIATED(boundary_conditions)) THEN
3643  CALL boundary_conditions_variable_get(solver_equations%BOUNDARY_CONDITIONS, &
3644  & field_variable,boundary_conditions_variable,err,error,*999)
3645  IF(ASSOCIATED(boundary_conditions_variable)) THEN
3646  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
3647  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE== &
3648  & field_node_based_interpolation) THEN
3649  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
3650  IF(ASSOCIATED(domain)) THEN
3651  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
3652  domain_nodes=>domain%TOPOLOGY%NODES
3653  IF(ASSOCIATED(domain_nodes)) THEN
3654  !Loop over the local nodes excluding the ghosts.
3655  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
3656 !!TODO \todo We should interpolate the geometric field here and the node position.
3657  DO dim_idx=1,number_of_dimensions
3658  !Default to version 1 of each node derivative
3659  local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP% &
3660  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(1)%VERSIONS(1)
3661  x(dim_idx)=geometric_parameters(local_ny)
3662  ENDDO !dim_idx
3663  !Loop over the derivatives
3664  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
3665  analytic_function_type=equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE
3666  global_deriv_index=domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)% &
3667  & global_derivative_index
3668  CALL diffusion_analyticfunctionsevaluate(equations_set, &
3669  & analytic_function_type,x,tangents,normal,current_time,variable_type, &
3670  & global_deriv_index,component_idx,analytic_parameters,materials_parameters, &
3671  & VALUE,err,error,*999)
3672  !Default to version 1 of each node derivative
3673  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
3674  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
3675  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
3676  & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
3677  boundary_condition_check_variable=boundary_conditions_variable% &
3678  & condition_types(local_ny)
3679  IF(boundary_condition_check_variable==boundary_condition_fixed) THEN
3680  CALL field_parameter_set_update_local_dof(dependent_field, &
3681  & variable_type,field_values_set_type,local_ny, &
3682  & VALUE,err,error,*999)
3683  ENDIF
3684 
3685  ! IF(variable_type==FIELD_U_VARIABLE_TYPE) THEN
3686  ! IF(DOMAIN_NODES%NODES(node_idx)%BOUNDARY_NODE) THEN
3687  !If we are a boundary node then set the analytic value on the boundary
3688  ! CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
3689  ! & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
3690  ! ENDIF
3691  ! ENDIF
3692  ENDDO !deriv_idx
3693  ENDDO !node_idx
3694  ELSE
3695  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
3696  ENDIF
3697  ELSE
3698  CALL flagerror("Domain topology is not associated.",err,error,*999)
3699  ENDIF
3700  ELSE
3701  CALL flagerror("Domain is not associated.",err,error,*999)
3702  ENDIF
3703  ELSE
3704  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
3705  ENDIF
3706  ENDDO !component_idx
3707  ELSE
3708  CALL flagerror("Boundary conditions variable is not associated",err,error,*999)
3709  ENDIF
3710  ELSE
3711  CALL flagerror("Boundary conditions are not associated",err,error,*999)
3712  ENDIF
3713  CALL field_parameter_set_update_start(dependent_field,variable_type, &
3714  & field_analytic_values_set_type,err,error,*999)
3715  CALL field_parameter_set_update_finish(dependent_field,variable_type, &
3716  & field_analytic_values_set_type,err,error,*999)
3717  CALL field_parameter_set_update_start(dependent_field,variable_type, &
3718  & field_values_set_type,err,error,*999)
3719  CALL field_parameter_set_update_finish(dependent_field,variable_type, &
3720  & field_values_set_type,err,error,*999)
3721  ELSE
3722  CALL flagerror("Field variable is not associated.",err,error,*999)
3723  ENDIF
3724 
3725  ! ENDDO !variable_idx
3726  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,&
3727  & field_values_set_type,geometric_parameters,err,error,*999)
3728  ELSE
3729  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
3730  ENDIF
3731  ELSE
3732  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
3733  ENDIF
3734  ELSE
3735  !CALL FlagError("Equations set analytic is not associated.",ERR,ERROR,*999)
3736  ENDIF
3737  ELSE
3738  CALL flagerror("Equations set is not associated.",err,error,*999)
3739  ENDIF
3740  ELSE
3741  CALL flagerror("Equations are not associated.",err,error,*999)
3742  END IF
3743  ! ELSE
3744  ! CALL FlagError("Solver equations are not associated.",ERR,ERROR,*999)
3745  ! END IF
3746  CALL field_parameter_set_update_start(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
3747  & field_values_set_type,err,error,*999)
3748  CALL field_parameter_set_update_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
3749  & field_values_set_type,err,error,*999)
3750  IF(control_loop%PROBLEM%SPECIFICATION(3)==problem_linear_source_diffusion_subtype)THEN
3752  IF(ASSOCIATED(equations_set)) THEN
3753  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
3754  equations_source=>equations_set%SOURCE
3755  IF(ASSOCIATED(equations_source)) THEN
3756  source_field=>equations_source%SOURCE_FIELD
3757  IF(ASSOCIATED(source_field)) THEN
3758  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
3759  IF(ASSOCIATED(geometric_field)) THEN
3760  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions, &
3761  & err,error,*999)
3762  NULLIFY(geometric_variable)
3763  CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
3764  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type, &
3765  & geometric_parameters,err,error,*999)
3766  variable_type=field_u_variable_type
3767  field_variable=>source_field%VARIABLE_TYPE_MAP(variable_type)%PTR
3768  IF(ASSOCIATED(field_variable)) THEN
3769  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
3770  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
3771  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
3772  IF(ASSOCIATED(domain)) THEN
3773  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
3774  domain_nodes=>domain%TOPOLOGY%NODES
3775  IF(ASSOCIATED(domain_nodes)) THEN
3776  !Loop over the local nodes excluding the ghosts.
3777  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
3778 !!TODO \todo We should interpolate the geometric field here and the node position.
3779  DO dim_idx=1,number_of_dimensions
3780  !Default to version 1 of each node derivative
3781  local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
3782  & nodes(node_idx)%DERIVATIVES(1)%VERSIONS(1)
3783  x(dim_idx)=geometric_parameters(local_ny)
3784  ENDDO !dim_idx
3785  !Loop over the derivatives
3786  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
3787  SELECT CASE(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE)
3789  value_source=-1*a1*exp(-1*current_time)*(x(1)*x(1)+x(2)*x(2)+x(3)*x(3)+6)
3790  CASE DEFAULT
3791  local_error="The analytic function type of "// &
3792  & trim(number_to_vstring(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE,"*", &
3793  & err,error))//" is invalid."
3794  CALL flagerror(local_error,err,error,*999)
3795  END SELECT
3796  !Default to version 1 of each node derivative
3797  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
3798  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
3799  CALL field_parameter_set_update_local_dof(source_field,field_u_variable_type, &
3800  & field_values_set_type,local_ny,value_source,err,error,*999)
3801  ENDDO !deriv_idx
3802  ENDDO !node_idx
3803  ELSE
3804  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
3805  ENDIF
3806  ELSE
3807  CALL flagerror("Domain topology is not associated.",err,error,*999)
3808  ENDIF
3809  ELSE
3810  CALL flagerror("Domain is not associated.",err,error,*999)
3811  ENDIF
3812  ELSE
3813  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
3814  ENDIF
3815  ENDDO !component_idx
3816  CALL field_parameter_set_update_start(source_field,field_u_variable_type,field_values_set_type, &
3817  & err,error,*999)
3818  CALL field_parameter_set_update_finish(source_field,field_u_variable_type,field_values_set_type, &
3819  & err,error,*999)
3820  ELSE
3821  CALL flagerror("Field variable is not associated.",err,error,*999)
3822  ENDIF
3823  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
3824  & geometric_parameters,err,error,*999)
3825  ELSE
3826  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
3827  ENDIF
3828  ELSE
3829  CALL flagerror("Equations set source field is not associated.",err,error,*999)
3830  ENDIF
3831  ENDIF
3832  ELSE
3833  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
3834  ENDIF
3835  ELSE
3836  CALL flagerror("Equations set is not associated.",err,error,*999)
3837  ENDIF
3838  ENDIF
3839  ENDDO !eqnset_idx
3840  ELSE
3841  CALL flagerror("Solver equations are not associated.",err,error,*999)
3842  END IF
3843  CASE DEFAULT
3844  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
3845  & " is not valid for a diffusion equation type of a classical field problem class."
3846  CALL flagerror(local_error,err,error,*999)
3847  END SELECT
3848  ELSE
3849  CALL flagerror("Problem is not associated.",err,error,*999)
3850  ENDIF
3851  ELSE
3852  CALL flagerror("Solver is not associated.",err,error,*999)
3853  ENDIF
3854  ELSE
3855  CALL flagerror("Control loop is not associated.",err,error,*999)
3856  ENDIF
3857 
3858  exits("Diffusion_PreSolveUpdateAnalyticValues")
3859  RETURN
3860 999 errors("Diffusion_PreSolveUpdateAnalyticValues",err,error)
3861  exits("Diffusion_PreSolveUpdateAnalyticValues")
3862  RETURN 1
3863 
3865 
3866  !
3867  !================================================================================================================================
3868  !
3870  SUBROUTINE diffusion_presolvealeupdatemesh(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
3872  !Argument variables
3873  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
3874  TYPE(solver_type), POINTER :: SOLVER
3875  INTEGER(INTG), INTENT(OUT) :: ERR
3876  TYPE(varying_string), INTENT(OUT) :: ERROR
3877  !Local Variables
3878  TYPE(field_type), POINTER :: GEOMETRIC_FIELD
3879  TYPE(solver_type), POINTER :: SOLVER_ALE_DIFFUSION
3880  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
3881  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
3882  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
3883  TYPE(varying_string) :: LOCAL_ERROR
3884 
3885  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT,ALPHA
3886  REAL(DP), POINTER :: MESH_DISPLACEMENT_VALUES(:)
3887 
3888  INTEGER(INTG) :: dof_number,TOTAL_NUMBER_OF_DOFS,NDOFS_TO_PRINT
3889 
3890  INTEGER(INTG) :: NUMBER_OF_DIMENSIONS
3891  INTEGER(INTG) :: INPUT_TYPE,INPUT_OPTION
3892  REAL(DP), POINTER :: INPUT_DATA1(:)
3893 
3894  enters("Diffusion_PreSolveALEUpdateMesh",err,error,*999)
3895 
3896  NULLIFY(solver_ale_diffusion)
3897  NULLIFY(solver_equations)
3898  NULLIFY(solver_mapping)
3899  NULLIFY(equations_set)
3900 
3901  IF(ASSOCIATED(control_loop)) THEN
3902  IF(control_loop%LOOP_TYPE==problem_control_time_loop_type) THEN
3903  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
3904  ELSE IF(control_loop%CONTROL_LOOP_LEVEL>1) THEN
3905  CALL control_loop_current_times_get(control_loop%PARENT_LOOP,current_time,time_increment,err,error,*999)
3906  ENDIF
3907  IF(ASSOCIATED(solver)) THEN
3908  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
3909  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
3910  CALL flagerror("Problem specification is not allocated.",err,error,*999)
3911  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
3912  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
3913  END IF
3914  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
3917  ! do nothing ???
3920  solver_equations=>solver%SOLVER_EQUATIONS
3921  IF(ASSOCIATED(solver_equations)) THEN
3922  solver_mapping=>solver_equations%SOLVER_MAPPING
3923  IF(ASSOCIATED(solver_mapping)) THEN
3924  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
3925  IF(ASSOCIATED(equations_set)) THEN
3926  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
3927  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
3928  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
3929  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
3930  & err,error,*999)
3931  END IF
3932  SELECT CASE(equations_set%SPECIFICATION(3))
3936  ! do nothing ???
3941  CALL write_string(general_output_type,"Diffusion update mesh ... ",err,error,*999)
3942  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
3943  IF(ASSOCIATED(geometric_field)) THEN
3944  !--- First, read mesh displacement values from file
3945 
3946  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
3947  & number_of_dimensions,err,error,*999)
3948 
3949  input_type=42
3950  input_option=2
3951  NULLIFY(input_data1)
3952  !CALL FIELD_PARAMETER_SET_DATA_GET(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, &
3953  !& FIELD_VALUES_SET_TYPE,INPUT_DATA1,ERR,ERROR,*999)
3954 ! CALL FLUID_MECHANICS_IO_READ_DATA(SOLVER_LINEAR_TYPE,INPUT_DATA1, &
3955 ! & NUMBER_OF_DIMENSIONS,INPUT_TYPE,INPUT_OPTION,CURRENT_TIME)
3956 
3957  NULLIFY(mesh_displacement_values)
3958  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type, &
3959  & field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
3960  IF(diagnostics1) THEN
3961  ndofs_to_print = SIZE(mesh_displacement_values,1)
3962  CALL write_string_vector(diagnostic_output_type,1,1,ndofs_to_print,ndofs_to_print,ndofs_to_print,&
3963  & mesh_displacement_values,'(" MESH_DISPLACEMENT_VALUES = ",3(X,E13.6))','3(3(X,E13.6))', &
3964  & err,error,*999)
3965  ENDIF
3966 
3967 ! CALL FLUID_MECHANICS_IO_READ_DATA(SOLVER_LINEAR_TYPE,INPUT_DATA1, &
3968 ! & NUMBER_OF_DIMENSIONS,INPUT_TYPE,INPUT_OPTION,CURRENT_TIME)
3969 
3970  total_number_of_dofs = geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR% &
3971  & total_number_of_dofs
3972 
3973  !--- Second, update geometric field
3974  DO dof_number=1,total_number_of_dofs
3975  CALL field_parameter_set_add_local_dof(geometric_field, &
3976  & field_u_variable_type,field_values_set_type,dof_number, &
3977  & mesh_displacement_values(dof_number), &
3978  & err,error,*999)
3979  END DO
3980  CALL field_parameter_set_update_start(geometric_field, &
3981  & field_u_variable_type, field_values_set_type,err,error,*999)
3982  CALL field_parameter_set_update_finish(geometric_field, &
3983  & field_u_variable_type, field_values_set_type,err,error,*999)
3984 
3985  !--- Third, use displacement values to calculate velocity values
3986  alpha=1.0_dp/time_increment
3987  CALL field_parameter_sets_copy(geometric_field,field_u_variable_type, &
3988  & field_mesh_displacement_set_type,field_mesh_velocity_set_type,alpha,err,error,*999)
3989  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type, &
3990  & field_mesh_displacement_set_type,mesh_displacement_values,err,error,*999)
3991  ELSE
3992  CALL flagerror("Geometric field is not associated.",err,error,*999)
3993  ENDIF
3994  CASE DEFAULT
3995  local_error="Equations set subtype " &
3996  & //trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
3997  & " is not valid for a diffusion equation type of a classical field problem class."
3998  CALL flagerror(local_error,err,error,*999)
3999  END SELECT
4000  ELSE
4001  CALL flagerror("Equations set is not associated.",err,error,*999)
4002  ENDIF
4003  ELSE
4004  CALL flagerror("Solver mapping is not associated.",err,error,*999)
4005  ENDIF
4006  ELSE
4007  CALL flagerror("Solver equations is not associated.",err,error,*999)
4008  ENDIF
4009  CASE DEFAULT
4010  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4011  & " is not valid for a diffusion equation type of a classical field problem class."
4012  CALL flagerror(local_error,err,error,*999)
4013  END SELECT
4014  ELSE
4015  CALL flagerror("Problem is not associated.",err,error,*999)
4016  ENDIF
4017  ELSE
4018  CALL flagerror("Solver is not associated.",err,error,*999)
4019  ENDIF
4020  ELSE
4021  CALL flagerror("Control loop is not associated.",err,error,*999)
4022  ENDIF
4023 
4024  exits("Diffusion_PreSolveALEUpdateMesh")
4025  RETURN
4026 999 errorsexits("Diffusion_PreSolveALEUpdateMesh",err,error)
4027  RETURN 1
4028 
4029  END SUBROUTINE diffusion_presolvealeupdatemesh
4030  !
4031  !================================================================================================================================
4032  !
4033  SUBROUTINE diffusion_presolvestorecurrentsolution(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4035  !Argument variables
4036  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4037  TYPE(solver_type), POINTER :: SOLVER
4038  INTEGER(INTG), INTENT(OUT) :: ERR
4039  TYPE(varying_string), INTENT(OUT) :: ERROR
4040 
4041  !Local Variables
4042  TYPE(solver_type), POINTER :: SOLVER_DIFFUSION_ONE
4043  TYPE(field_type), POINTER :: DEPENDENT_FIELD_DIFFUSION_ONE
4044  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS_DIFFUSION_ONE
4045  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING_DIFFUSION_ONE
4046  TYPE(equations_set_type), POINTER :: EQUATIONS_SET_DIFFUSION_ONE
4047  TYPE(varying_string) :: LOCAL_ERROR
4048 
4049  INTEGER(INTG) :: NUMBER_OF_COMPONENTS_DEPENDENT_FIELD_DIFFUSION_ONE
4050  INTEGER(INTG) :: I
4051 
4052  enters("Diffusion_PreSolveStoreCurrentSolution",err,error,*999)
4053 
4054  IF(ASSOCIATED(control_loop)) THEN
4055 
4056  NULLIFY(solver_diffusion_one)
4057 
4058  IF(ASSOCIATED(solver)) THEN
4059  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4060  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4061  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4062  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4063  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
4064  END IF
4065  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4067  ! do nothing ???
4069  ! do nothing ???
4071  ! do nothing ???
4073  IF(solver%GLOBAL_NUMBER==1) THEN
4074  !--- Get the dependent field of the diffusion-one equations
4075  CALL write_string(general_output_type,"Store value diffusion-one dependent field at time, t ... ",err,error,*999)
4076  CALL solvers_solver_get(solver%SOLVERS,1,solver_diffusion_one,err,error,*999)
4077  solver_equations_diffusion_one=>solver_diffusion_one%SOLVER_EQUATIONS
4078  IF(ASSOCIATED(solver_equations_diffusion_one)) THEN
4079  solver_mapping_diffusion_one=>solver_equations_diffusion_one%SOLVER_MAPPING
4080  IF(ASSOCIATED(solver_mapping_diffusion_one)) THEN
4081  equations_set_diffusion_one=>solver_mapping_diffusion_one%EQUATIONS_SETS(1)%PTR
4082  IF(ASSOCIATED(equations_set_diffusion_one)) THEN
4083  dependent_field_diffusion_one=>equations_set_diffusion_one%DEPENDENT%DEPENDENT_FIELD
4084  IF(ASSOCIATED(dependent_field_diffusion_one)) THEN
4085  CALL field_number_of_components_get(dependent_field_diffusion_one, &
4086  & field_u_variable_type,number_of_components_dependent_field_diffusion_one,err,error,*999)
4087  ELSE
4088  CALL flagerror("DEPENDENT_FIELD_DIFFUSION_ONE is not associated.",err,error,*999)
4089  END IF
4090  ELSE
4091  CALL flagerror("Diffusion-one equations set is not associated.",err,error,*999)
4092  END IF
4093  ELSE
4094  CALL flagerror("Diffusion-one solver mapping is not associated.",err,error,*999)
4095  END IF
4096  ELSE
4097  CALL flagerror("Diffusion-one solver equations are not associated.",err,error,*999)
4098  END IF
4099 
4100  !--- Copy the current time value parameters set from diffusion-one's dependent field
4101  DO i=1,number_of_components_dependent_field_diffusion_one
4102  CALL field_parameterstofieldparameterscopy(dependent_field_diffusion_one, &
4103  & field_u_variable_type,field_values_set_type,i,dependent_field_diffusion_one, &
4104  & field_u_variable_type,field_previous_values_set_type,i,err,error,*999)
4105  END DO
4106 
4107 
4108 ! IF(DIAGNOSTICS3) THEN
4109 ! NULLIFY( DUMMY_VALUES2 )
4110 ! CALL FIELD_PARAMETER_SET_DATA_GET(DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE, &
4111 ! & FIELD_VALUES_SET_TYPE,DUMMY_VALUES2,ERR,ERROR,*999)
4112 ! NDOFS_TO_PRINT = SIZE(DUMMY_VALUES2,1)
4113 ! CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,NDOFS_TO_PRINT,NDOFS_TO_PRINT,NDOFS_TO_PRINT,DUMMY_VALUES2, &
4114 ! & '(" DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE = ",4(X,E13.6))',&
4115 ! & '4(4(X,E13.6))',ERR,ERROR,*999)
4116 ! CALL FIELD_PARAMETER_SET_DATA_RESTORE(DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE, &
4117 ! & FIELD_VALUES_SET_TYPE,DUMMY_VALUES2,ERR,ERROR,*999)
4118 ! ENDIF
4119 
4120  END IF
4122  IF(solver%GLOBAL_NUMBER==2) THEN
4123  !--- Get the dependent field of the diffusion equations
4124  CALL write_string(general_output_type,"Store value of diffusion solution &
4125  & (dependent field - V variable_type) at time, t ... ",err,error,*999)
4126  CALL solvers_solver_get(solver%SOLVERS,2,solver_diffusion_one,err,error,*999)
4127  solver_equations_diffusion_one=>solver_diffusion_one%SOLVER_EQUATIONS
4128  IF(ASSOCIATED(solver_equations_diffusion_one)) THEN
4129  solver_mapping_diffusion_one=>solver_equations_diffusion_one%SOLVER_MAPPING
4130  IF(ASSOCIATED(solver_mapping_diffusion_one)) THEN
4131  equations_set_diffusion_one=>solver_mapping_diffusion_one%EQUATIONS_SETS(1)%PTR
4132  IF(ASSOCIATED(equations_set_diffusion_one)) THEN
4133  dependent_field_diffusion_one=>equations_set_diffusion_one%DEPENDENT%DEPENDENT_FIELD
4134  IF(ASSOCIATED(dependent_field_diffusion_one)) THEN
4135  CALL field_number_of_components_get(dependent_field_diffusion_one, &
4136  & field_v_variable_type,number_of_components_dependent_field_diffusion_one,err,error,*999)
4137  ELSE
4138  CALL flagerror("DEPENDENT_FIELD_DIFFUSION_ONE is not associated.",err,error,*999)
4139  END IF
4140  ELSE
4141  CALL flagerror("Diffusion equations set is not associated.",err,error,*999)
4142  END IF
4143  ELSE
4144  CALL flagerror("Diffusion solver mapping is not associated.",err,error,*999)
4145  END IF
4146  ELSE
4147  CALL flagerror("Diffusion solver equations are not associated.",err,error,*999)
4148  END IF
4149 
4150  !--- Copy the current time value parameters set from diffusion-one's dependent field
4151  DO i=1,number_of_components_dependent_field_diffusion_one
4152  CALL field_parameterstofieldparameterscopy(dependent_field_diffusion_one, &
4153  & field_v_variable_type,field_values_set_type,i,dependent_field_diffusion_one, &
4154  & field_v_variable_type,field_previous_values_set_type,i,err,error,*999)
4155  END DO
4156 
4157 
4158 ! IF(DIAGNOSTICS3) THEN
4159 ! NULLIFY( DUMMY_VALUES2 )
4160 ! CALL FIELD_PARAMETER_SET_DATA_GET(DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE, &
4161 ! & FIELD_VALUES_SET_TYPE,DUMMY_VALUES2,ERR,ERROR,*999)
4162 ! NDOFS_TO_PRINT = SIZE(DUMMY_VALUES2,1)
4163 ! CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,NDOFS_TO_PRINT,NDOFS_TO_PRINT,NDOFS_TO_PRINT,DUMMY_VALUES2, &
4164 ! & '(" DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE = ",4(X,E13.6))',&
4165 ! & '4(4(X,E13.6))',ERR,ERROR,*999)
4166 ! CALL FIELD_PARAMETER_SET_DATA_RESTORE(DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE, &
4167 ! & FIELD_VALUES_SET_TYPE,DUMMY_VALUES2,ERR,ERROR,*999)
4168 ! ENDIF
4169 
4170  END IF
4171  CASE DEFAULT
4172  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4173  & " is not valid for a diffusion equation type of a classical field problem class."
4174  CALL flagerror(local_error,err,error,*999)
4175  END SELECT
4176  ELSE
4177  CALL flagerror("Problem is not associated.",err,error,*999)
4178  ENDIF
4179  ELSE
4180  CALL flagerror("Solver is not associated.",err,error,*999)
4181  ENDIF
4182  ELSE
4183  CALL flagerror("Control loop is not associated.",err,error,*999)
4184  ENDIF
4185 
4186  exits("Diffusion_PreSolveStoreCurrentSolution")
4187  RETURN
4188 999 errors("Diffusion_PreSolveStoreCurrentSolution",err,error)
4189  exits("Diffusion_PreSolveStoreCurrentSolution")
4190  RETURN 1
4191 
4193 
4194  !
4195  !================================================================================================================================
4196  !
4197  SUBROUTINE diffusion_presolvegetsourcevalue(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4199  !Argument variables
4200  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4201  TYPE(solver_type), POINTER :: SOLVER
4202  INTEGER(INTG), INTENT(OUT) :: ERR
4203  TYPE(varying_string), INTENT(OUT) :: ERROR
4204 
4205  !Local Variables
4206  TYPE(solver_type), POINTER :: SOLVER_DIFFUSION_ONE, SOLVER_DIFFUSION_TWO
4207  TYPE(field_type), POINTER :: DEPENDENT_FIELD_DIFFUSION_TWO, SOURCE_FIELD_DIFFUSION_ONE
4208  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS_DIFFUSION_ONE, SOLVER_EQUATIONS_DIFFUSION_TWO
4209  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING_DIFFUSION_ONE, SOLVER_MAPPING_DIFFUSION_TWO
4210  TYPE(equations_set_type), POINTER :: EQUATIONS_SET_DIFFUSION_ONE, EQUATIONS_SET_DIFFUSION_TWO
4211  TYPE(varying_string) :: LOCAL_ERROR
4212 
4213  INTEGER(INTG) :: NUMBER_OF_COMPONENTS_DEPENDENT_FIELD_DIFFUSION_TWO,NUMBER_OF_COMPONENTS_SOURCE_FIELD_DIFFUSION_ONE
4214  INTEGER(INTG) :: I
4215 
4216 
4217  enters("Diffusion_PreSolveGetSourceValue",err,error,*999)
4218 
4219  IF(ASSOCIATED(control_loop)) THEN
4220 
4221  NULLIFY(solver_diffusion_one)
4222  NULLIFY(solver_diffusion_two)
4223 
4224  IF(ASSOCIATED(solver)) THEN
4225  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4226  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4227  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4228  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4229  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
4230  END IF
4231  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4233  ! do nothing ???
4235  ! do nothing ???
4237  ! do nothing ???
4239  IF(solver%GLOBAL_NUMBER==1) THEN
4240  !--- Get the dependent field of the diffusion_two equations
4241  CALL write_string(general_output_type,"Update diffusion-one source field ... ",err,error,*999)
4242  CALL solvers_solver_get(solver%SOLVERS,2,solver_diffusion_two,err,error,*999)
4243  solver_equations_diffusion_two=>solver_diffusion_two%SOLVER_EQUATIONS
4244  IF(ASSOCIATED(solver_equations_diffusion_two)) THEN
4245  solver_mapping_diffusion_two=>solver_equations_diffusion_two%SOLVER_MAPPING
4246  IF(ASSOCIATED(solver_mapping_diffusion_two)) THEN
4247  equations_set_diffusion_two=>solver_mapping_diffusion_two%EQUATIONS_SETS(1)%PTR
4248  IF(ASSOCIATED(equations_set_diffusion_two)) THEN
4249  dependent_field_diffusion_two=>equations_set_diffusion_two%DEPENDENT%DEPENDENT_FIELD
4250  IF(ASSOCIATED(dependent_field_diffusion_two)) THEN
4251  CALL field_number_of_components_get(dependent_field_diffusion_two, &
4252  & field_u_variable_type,number_of_components_dependent_field_diffusion_two,err,error,*999)
4253  ELSE
4254  CALL flagerror("DEPENDENT_FIELD_DIFFUSION_TWO is not associated.",err,error,*999)
4255  END IF
4256  ELSE
4257  CALL flagerror("Diffusion-two equations set is not associated.",err,error,*999)
4258  END IF
4259  ELSE
4260  CALL flagerror("Diffusion-two solver mapping is not associated.",err,error,*999)
4261  END IF
4262  ELSE
4263  CALL flagerror("Diffusion-two solver equations are not associated.",err,error,*999)
4264  END IF
4265 
4266 
4267  !--- Get the source field for the diffusion_one equations
4268  CALL solvers_solver_get(solver%SOLVERS,1,solver_diffusion_one,err,error,*999)
4269  solver_equations_diffusion_one=>solver_diffusion_one%SOLVER_EQUATIONS
4270  IF(ASSOCIATED(solver_equations_diffusion_one)) THEN
4271  solver_mapping_diffusion_one=>solver_equations_diffusion_one%SOLVER_MAPPING
4272  IF(ASSOCIATED(solver_mapping_diffusion_one)) THEN
4273  equations_set_diffusion_one=>solver_mapping_diffusion_one%EQUATIONS_SETS(1)%PTR
4274  IF(ASSOCIATED(equations_set_diffusion_one)) THEN
4275  source_field_diffusion_one=>equations_set_diffusion_one%SOURCE%SOURCE_FIELD
4276  IF(ASSOCIATED(source_field_diffusion_one)) THEN
4277  CALL field_number_of_components_get(source_field_diffusion_one, &
4278  & field_u_variable_type,number_of_components_source_field_diffusion_one,err,error,*999)
4279  ELSE
4280  CALL flagerror("SOURCE_FIELD_DIFFUSION_ONE is not associated.",err,error,*999)
4281  END IF
4282  ELSE
4283  CALL flagerror("Diffusion-one equations set is not associated.",err,error,*999)
4284  END IF
4285  ELSE
4286  CALL flagerror("Diffusion-one solver mapping is not associated.",err,error,*999)
4287  END IF
4288  ELSE
4289  CALL flagerror("Diffusion-one solver equations are not associated.",err,error,*999)
4290  END IF
4291 
4292  !--- Copy the result from diffusion-two's dependent field to diffusion-one's source field
4293  IF(number_of_components_source_field_diffusion_one==number_of_components_dependent_field_diffusion_two) THEN
4294  DO i=1,number_of_components_source_field_diffusion_one
4295  CALL field_parameterstofieldparameterscopy(dependent_field_diffusion_two, &
4296  & field_u_variable_type,field_values_set_type,i,source_field_diffusion_one, &
4297  & field_u_variable_type,field_values_set_type,i,err,error,*999)
4298 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!FIELDMESHDISPLACEMENTTYPE needs to be changed to appropriate type for this problem
4299  END DO
4300  ELSE
4301  local_error="Number of components of diffusion-two dependent field "// &
4302  & "is not consistent with diffusion-one-equation source field."
4303  CALL flagerror(local_error,err,error,*999)
4304  END IF
4305 
4306 ! IF(DIAGNOSTICS3) THEN
4307 ! NULLIFY( DUMMY_VALUES2 )
4308 ! CALL FIELD_PARAMETER_SET_DATA_GET(DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE, &
4309 ! & FIELD_VALUES_SET_TYPE,DUMMY_VALUES2,ERR,ERROR,*999)
4310 ! NDOFS_TO_PRINT = SIZE(DUMMY_VALUES2,1)
4311 ! CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,NDOFS_TO_PRINT,NDOFS_TO_PRINT,NDOFS_TO_PRINT,DUMMY_VALUES2, &
4312 ! & '(" DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE = ",4(X,E13.6))',&
4313 ! & '4(4(X,E13.6))',ERR,ERROR,*999)
4314 ! CALL FIELD_PARAMETER_SET_DATA_RESTORE(DEPENDENT_FIELD_FINITE_ELASTICITY,FIELD_U_VARIABLE_TYPE, &
4315 ! & FIELD_VALUES_SET_TYPE,DUMMY_VALUES2,ERR,ERROR,*999)
4316 ! ENDIF
4317 
4318  END IF
4319  CASE DEFAULT
4320  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4321  & " is not valid for a diffusion equation type of a classical field problem class."
4322  CALL flagerror(local_error,err,error,*999)
4323  END SELECT
4324  ELSE
4325  CALL flagerror("Problem is not associated.",err,error,*999)
4326  ENDIF
4327  ELSE
4328  CALL flagerror("Solver is not associated.",err,error,*999)
4329  ENDIF
4330  ELSE
4331  CALL flagerror("Control loop is not associated.",err,error,*999)
4332  ENDIF
4333 
4334  exits("Diffusion_PreSolveGetSourceValue")
4335  RETURN
4336 999 errorsexits("Diffusion_PreSolveGetSourceValue",err,error)
4337  RETURN 1
4338 
4339  END SUBROUTINE diffusion_presolvegetsourcevalue
4340  !
4341  !================================================================================================================================
4342  !
4344  SUBROUTINE diffusion_equation_post_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4346  !Argument variables
4347  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4348  TYPE(solver_type), POINTER :: SOLVER
4349  INTEGER(INTG), INTENT(OUT) :: ERR
4350  TYPE(varying_string), INTENT(OUT) :: ERROR
4351  !Local Variables
4352  TYPE(solver_type), POINTER :: SOLVER2
4353  TYPE(varying_string) :: LOCAL_ERROR
4354 
4355  enters("DIFFUSION_EQUATION_POST_SOLVE",err,error,*999)
4356 
4357  NULLIFY(solver2)
4358 
4359  IF(ASSOCIATED(control_loop)) THEN
4360  IF(ASSOCIATED(solver)) THEN
4361  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4362  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4363  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4364  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4365  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
4366  END IF
4367  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4370  !CALL DIFFUSION_EQUATION_POST_SOLVE_OUTPUT_DATA(CONTROL_LOOP,SOLVER,ERR,ERROR,*999)
4372  ! do nothing ???
4374  CALL flagerror("Not implemented.",err,error,*999)
4375  CASE DEFAULT
4376  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4377  & " is not valid for a diffusion type of a classical field problem class."
4378  CALL flagerror(local_error,err,error,*999)
4379  END SELECT
4380  ELSE
4381  CALL flagerror("Problem is not associated.",err,error,*999)
4382  ENDIF
4383  ELSE
4384  CALL flagerror("Solver is not associated.",err,error,*999)
4385  ENDIF
4386  ELSE
4387  CALL flagerror("Control loop is not associated.",err,error,*999)
4388  ENDIF
4389 
4390  exits("DIFFUSION_EQUATION_POST_SOLVE")
4391  RETURN
4392 999 errorsexits("DIFFUSION_EQUATION_POST_SOLVE",err,error)
4393  RETURN 1
4394 
4395  END SUBROUTINE diffusion_equation_post_solve
4396 
4397  !
4398  !================================================================================================================================
4399  !
4400 
4402  SUBROUTINE diffusion_equation_post_solve_output_data(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
4404  !Argument variables
4405  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
4406  TYPE(solver_type), POINTER :: SOLVER
4407  INTEGER(INTG), INTENT(OUT) :: ERR
4408  TYPE(varying_string), INTENT(OUT) :: ERROR
4409  !Local Variables
4410  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
4411  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
4412  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4413  TYPE(varying_string) :: LOCAL_ERROR
4414 
4415  REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
4416  INTEGER(INTG) :: EQUATIONS_SET_IDX,CURRENT_LOOP_ITERATION,OUTPUT_ITERATION_NUMBER
4417 
4418  CHARACTER(14) :: FILE
4419  CHARACTER(14) :: OUTPUT_FILE
4420 
4421  enters("DIFFUSION_EQUATION_POST_SOLVE_OUTPUT_DATA",err,error,*999)
4422 
4423  IF(ASSOCIATED(control_loop)) THEN
4424 ! write(*,*)'CURRENT_TIME = ',CURRENT_TIME
4425 ! write(*,*)'TIME_INCREMENT = ',TIME_INCREMENT
4426  IF(ASSOCIATED(solver)) THEN
4427  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
4428  IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
4429  CALL flagerror("Problem specification is not allocated.",err,error,*999)
4430  ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
4431  CALL flagerror("Problem specification must have three entries for a Diffusion problem.",err,error,*999)
4432  END IF
4433  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
4435  CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
4436  solver_equations=>solver%SOLVER_EQUATIONS
4437  IF(ASSOCIATED(solver_equations)) THEN
4438  solver_mapping=>solver_equations%SOLVER_MAPPING
4439  IF(ASSOCIATED(solver_mapping)) THEN
4440  !Make sure the equations sets are up to date
4441  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
4442  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
4443 
4444  current_loop_iteration=control_loop%TIME_LOOP%ITERATION_NUMBER
4445  output_iteration_number=control_loop%TIME_LOOP%OUTPUT_NUMBER
4446 
4447  IF(output_iteration_number/=0) THEN
4448  IF(control_loop%TIME_LOOP%CURRENT_TIME<=control_loop%TIME_LOOP%STOP_TIME) THEN
4449  IF(current_loop_iteration<10) THEN
4450  WRITE(output_file,'("TIME_STEP_000",I0)') current_loop_iteration
4451  ELSE IF(current_loop_iteration<100) THEN
4452  WRITE(output_file,'("TIME_STEP_00",I0)') current_loop_iteration
4453  ELSE IF(current_loop_iteration<1000) THEN
4454  WRITE(output_file,'("TIME_STEP_0",I0)') current_loop_iteration
4455  ELSE IF(current_loop_iteration<10000) THEN
4456  WRITE(output_file,'("TIME_STEP_",I0)') current_loop_iteration
4457  END IF
4458  file=output_file
4459  ! FILE="TRANSIENT_OUTPUT"
4460 !!!!!!!!ADAPT THIS TO WORK WITH DIFFUSION AND NOT JUST FLUID MECHANICS
4461 ! METHOD="FORTRAN"
4462 ! EXPORT_FIELD=.TRUE.
4463 ! IF(EXPORT_FIELD) THEN
4464 ! IF(MOD(CURRENT_LOOP_ITERATION,OUTPUT_ITERATION_NUMBER)==0) THEN
4465 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"...",ERR,ERROR,*999)
4466 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"Now export fields... ",ERR,ERROR,*999)
4467 ! CALL FLUID_MECHANICS_IO_WRITE_CMGUI(EQUATIONS_SET%REGION,EQUATIONS_SET%GLOBAL_NUMBER,FILE, &
4468 ! & ERR,ERROR,*999)
4469 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,OUTPUT_FILE,ERR,ERROR,*999)
4470 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,"...",ERR,ERROR,*999)
4471 ! ENDIF
4472 ! ENDIF
4473 
4474  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
4475  IF(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE==equations_set_diffusion_equation_two_dim_1 .OR. &
4476  & equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE== &
4478  CALL analyticanalysis_output(equations_set%DEPENDENT%DEPENDENT_FIELD,file,err,error,*999)
4479  ENDIF
4480  ENDIF
4481  ENDIF
4482  ENDIF
4483  ENDDO
4484  ENDIF
4485  ENDIF
4487  ! do nothing ???
4488  CALL flagerror("Not implemented.",err,error,*999)
4489  CASE DEFAULT
4490  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
4491  & " is not valid for a diffusion equation type of a classical field problem class."
4492  CALL flagerror(local_error,err,error,*999)
4493  END SELECT
4494  ELSE
4495  CALL flagerror("Problem is not associated.",err,error,*999)
4496  ENDIF
4497  ELSE
4498  CALL flagerror("Solver is not associated.",err,error,*999)
4499  ENDIF
4500  ELSE
4501  CALL flagerror("Control loop is not associated.",err,error,*999)
4502  ENDIF
4503  exits("DIFFUSION_EQUATION_POST_SOLVE_OUTPUT_DATA")
4504  RETURN
4505 999 errorsexits("DIFFUSION_EQUATION_POST_SOLVE_OUTPUT_DATA",err,error)
4506  RETURN 1
4508 
4509  !
4510  !================================================================================================================================
4511  !
4513  SUBROUTINE diffusion_equation_finite_element_calculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
4515  !Argument variables
4516  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
4517  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
4518  INTEGER(INTG), INTENT(OUT) :: ERR
4519  TYPE(varying_string), INTENT(OUT) :: ERROR
4520  !Local Variables
4521  INTEGER(INTG) FIELD_VAR_TYPE,mh,mhs,ms,ng,nh,nhs,ni,nj,ns,my_compartment,Ncompartments,imatrix,num_var_count
4522  INTEGER(INTG) :: MESH_COMPONENT_1, MESH_COMPONENT_2
4523  REAL(DP) :: C_PARAM,K_PARAM,RWG,SUM,PGMJ(3),PGNJ(3),A_PARAM,COUPLING_PARAM,PGM,PGN
4524  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
4525  TYPE(basis_type), POINTER :: DEPENDENT_BASIS_1, DEPENDENT_BASIS_2
4526  TYPE(equations_type), POINTER :: EQUATIONS
4527  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
4528  TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
4529  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
4530  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
4531  TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
4532  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
4533  TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
4534  TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
4535  TYPE(equations_matrix_type), POINTER :: DAMPING_MATRIX,STIFFNESS_MATRIX
4536  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD,SOURCE_FIELD,EQUATIONS_SET_FIELD
4537  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
4538  TYPE(field_variable_ptr_type) :: FIELD_VARIABLES(99)
4539  TYPE(equations_matrix_ptr_type) :: COUPLING_MATRICES(99)
4540  INTEGER(INTG) :: FIELD_VAR_TYPES(99)
4541  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
4542  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME_1, QUADRATURE_SCHEME_2
4543  TYPE(varying_string) :: LOCAL_ERROR
4544  TYPE(field_interpolation_parameters_type), POINTER :: ADVEC_DIFF_DEPENDENT_CURRENT_INTERPOLATION_PARAMETERS, &
4545  & ADVEC_DIFF_DEPENDENT_PREVIOUS_INTERPOLATION_PARAMETERS,DIFFUSION_DEPENDENT_PREVIOUS_INTERPOLATION_PARAMETERS
4546  TYPE(field_interpolated_point_type), POINTER :: ADVEC_DIFF_DEPENDENT_CURRENT_INTERPOLATED_POINT, &
4547  & ADVEC_DIFF_DEPENDENT_PREVIOUS_INTERPOLATED_POINT,DIFFUSION_DEPENDENT_PREVIOUS_INTERPOLATED_POINT
4548  INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:)
4549 
4550  enters("DIFFUSION_EQUATION_FINITE_ELEMENT_CALCULATE",err,error,*999)
4551 
4552  IF(ASSOCIATED(equations_set)) THEN
4553  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
4554  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
4555  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
4556  CALL flagerror("Equations set specification must have three entries for a diffusion type equations set.", &
4557  & err,error,*999)
4558  END IF
4559  equations=>equations_set%EQUATIONS
4560  IF(ASSOCIATED(equations)) THEN
4561  SELECT CASE(equations_set%SPECIFICATION(3))
4566  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
4567  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
4568  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
4569  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
4570  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
4571  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
4572  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
4573  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4574  source_field=>equations%INTERPOLATION%SOURCE_FIELD
4575  ENDIF
4576  equations_matrices=>equations%EQUATIONS_MATRICES
4577  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
4578  stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
4579  damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
4580  rhs_vector=>equations_matrices%RHS_VECTOR
4581  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
4582  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
4583  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
4584  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
4585  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4586  source_vector=>equations_matrices%SOURCE_VECTOR
4587  ENDIF
4588  equations_mapping=>equations%EQUATIONS_MAPPING
4589  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
4590  field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
4591  field_var_type=field_variable%VARIABLE_TYPE
4592  geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
4593  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4594  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4595  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4596  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4597  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
4598  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4599  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4600  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4601  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4602  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
4603  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
4604  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
4605  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
4606  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4607  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4608  & source_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4609  ENDIF
4610  IF(equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4611  advec_diff_dependent_current_interpolation_parameters=> &
4612  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_u_variable_type)%PTR
4613  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
4614  & advec_diff_dependent_current_interpolation_parameters,err,error,*999)
4615  advec_diff_dependent_current_interpolated_point=> &
4616  & equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_u_variable_type)%PTR
4617  diffusion_dependent_previous_interpolation_parameters=> &
4618  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_v_variable_type)%PTR
4619  CALL field_interpolation_parameters_element_get(field_previous_values_set_type,element_number, &
4620  & diffusion_dependent_previous_interpolation_parameters,err,error,*999)
4621  diffusion_dependent_previous_interpolated_point=> &
4622  & equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_v_variable_type)%PTR
4623  ENDIF
4624  !Loop over gauss points
4625  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
4626  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4627  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
4628  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
4629  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
4630  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4631  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
4632  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
4633  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
4634  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
4635  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
4636  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4637  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4638  & source_interp_point(field_u_variable_type)%PTR,err,error,*999)
4639  ENDIF
4640  !Calculate RWG.
4641 !!TODO: Think about symmetric problems.
4642  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
4643  & quadrature_scheme%GAUSS_WEIGHTS(ng)
4644  !Loop over field components
4645  mhs=0
4646  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
4647  !Loop over element rows
4648  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4649  mhs=mhs+1
4650  nhs=0
4651  IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX) THEN
4652  !Loop over element columns
4653  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
4654  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4655  nhs=nhs+1
4656  IF(stiffness_matrix%UPDATE_MATRIX) THEN
4657  sum=0.0_dp
4658  DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
4659  pgmj(nj)=0.0_dp
4660  pgnj(nj)=0.0_dp
4661  DO ni=1,dependent_basis%NUMBER_OF_XI
4662  pgmj(nj)=pgmj(nj)+ &
4663  & quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)* &
4664  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4665  pgnj(nj)=pgnj(nj)+ &
4666  & quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)* &
4667  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4668  ENDDO !ni
4669  k_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4670  & values(nj,no_part_deriv)
4671  sum=sum+k_param*pgmj(nj)*pgnj(nj)
4672  ENDDO !nj
4673  IF(equations_set%SPECIFICATION(3)==equations_set_no_source_diffusion_subtype .OR. &
4674  & equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
4675  & equations_set%SPECIFICATION(3)==equations_set_no_source_ale_diffusion_subtype .OR. &
4676  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
4677  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4678  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg
4679  ELSEIF(equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
4680  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype) THEN
4681  a_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4682  & values(geometric_variable%NUMBER_OF_COMPONENTS,no_part_deriv)
4683  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg- &
4684  & a_param*quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
4685  & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*rwg
4686  ENDIF
4687  ENDIF
4688  IF(damping_matrix%UPDATE_MATRIX) THEN
4689  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+ &
4690  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
4691  & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*rwg
4692  ENDIF
4693  ENDDO !ns
4694  ENDDO !nh
4695  ENDIF
4696  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
4697  ENDDO !ms
4698  ENDDO !mh
4699  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
4700  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
4701  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
4702  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype) THEN
4703  IF(source_vector%UPDATE_VECTOR) THEN
4704  c_param=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1, no_part_deriv)
4705  mhs=0
4706  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
4707  !Loop over element rows
4708  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4709  mhs=mhs+1
4710  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+ &
4711  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*c_param*rwg
4712  ENDDO !ms
4713  ENDDO !mh
4714  ENDIF
4715  ELSEIF(equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4716  IF(source_vector%UPDATE_VECTOR) THEN
4717  !The value of the source term is +0.5*(C_1^{t}+C_1_{t+1}-C_2^{t})
4718  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng, &
4719  & advec_diff_dependent_current_interpolated_point,err,error,*999)
4720  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng, &
4721  & diffusion_dependent_previous_interpolated_point,err,error,*999)
4722  write(*,*) advec_diff_dependent_current_interpolated_point%VALUES(1,no_part_deriv)
4723  write(*,*) diffusion_dependent_previous_interpolated_point%VALUES(1,no_part_deriv)
4724  c_param=0.5_dp*advec_diff_dependent_current_interpolated_point%VALUES(1,no_part_deriv)- &
4725  & diffusion_dependent_previous_interpolated_point%VALUES(1,no_part_deriv)
4726  ! C_PARAM_1_T= EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_POINT(FIELD_U_VARIABLE_TYPE)%PTR%VALUES(1, NO_PART_DERIV)!<This is the value of the solution from the advection-diffusion equation at time T
4727  ! C_PARAM_1_TPLUSONE= !<This is the value of the solution from the advection-diffusion equation at time T+deltaT
4728  ! C_PARAM_2_T= EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_POINT(FIELD_V_VARIABLE_TYPE)%PTR%VALUES(1, NO_PART_DERIV)!<This is the value of the solution from the diffusion equation at time T
4729  ! C_PARAM=C_PARAM_1_T+C_PARAM_1_TPLUSONE+C_PARAM_2_T
4730  mhs=0
4731  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
4732  !Loop over element rows
4733  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4734  mhs=mhs+1
4735  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+ &
4736  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*c_param*rwg
4737  ENDDO !ms
4738  ENDDO !mh
4739  ENDIF
4740  ENDIF
4741  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
4742  ENDDO !ng
4743 
4744  IF(equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4745  advec_diff_dependent_previous_interpolation_parameters=> &
4746  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_u_variable_type)%PTR
4747  CALL field_interpolation_parameters_element_get(field_previous_values_set_type,element_number, &
4748  & advec_diff_dependent_previous_interpolation_parameters,err,error,*999)
4749  advec_diff_dependent_previous_interpolated_point=> &
4750  & equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_u_variable_type)%PTR
4751  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
4752  IF(source_vector%UPDATE_VECTOR) THEN
4753  !The value of the source term is +0.5*(C_1^{t}+C_1_{t+1}-C_2^{t})
4754  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng, &
4755  & advec_diff_dependent_previous_interpolated_point,err,error,*999)
4756  write(*,*) advec_diff_dependent_previous_interpolated_point%VALUES(1,no_part_deriv)
4757  c_param=0.5_dp*advec_diff_dependent_previous_interpolated_point%VALUES(1,no_part_deriv)
4758  mhs=0
4759  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
4760  !Loop over element rows
4761  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4762  mhs=mhs+1
4763  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+ &
4764  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*c_param*rwg
4765  ENDDO !ms
4766  ENDDO !mh
4767  ENDIF
4768  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
4769  ENDDO !ng
4770  ENDIF
4771 
4772  !Scale factor adjustment
4773  IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
4774  CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
4775  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
4776  mhs=0
4777  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
4778  !Loop over element rows
4779  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4780  mhs=mhs+1
4781  nhs=0
4782  IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX) THEN
4783  !Loop over element columns
4784  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
4785  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4786  nhs=nhs+1
4787  IF(stiffness_matrix%UPDATE_MATRIX) THEN
4788  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
4789  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
4790  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
4791  ENDIF
4792  IF(damping_matrix%UPDATE_MATRIX) THEN
4793  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
4794  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
4795  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
4796  ENDIF
4797  ENDDO !ns
4798  ENDDO !nh
4799  ENDIF
4800  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
4801  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
4802  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_diffusion_subtype .OR. &
4803  & equations_set%SPECIFICATION(3)==equations_set_linear_source_diffusion_subtype .OR. &
4804  & equations_set%SPECIFICATION(3)==equations_set_constant_source_ale_diffusion_subtype .OR. &
4805  & equations_set%SPECIFICATION(3)==equations_set_linear_source_ale_diffusion_subtype .OR. &
4806  & equations_set%SPECIFICATION(3)==equations_set_coupled_source_diffusion_advec_diffusion_subtype) THEN
4807  IF(source_vector%UPDATE_VECTOR) source_vector%ELEMENT_VECTOR%VECTOR(mhs)= &
4808  & source_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
4809  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
4810  ENDIF
4811  ENDDO !ms
4812  ENDDO !mh
4813  ENDIF
4814 !!!!!!!!!!!!!!!MULTI-COMPARTMENT DIFFUSION - PROTOTYPE FOR OTHER MULTI-COMPARTMENT MODELS IN FUTURE
4815 !!!!!!!!!!!!!!!HAS BEEN SEPARATED HERE FOR EASE OF DEVELOPMENT & READABILITY OF THIS NEW FEATURE
4817 
4818  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
4819  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
4820  materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
4821  source_field=>equations%INTERPOLATION%SOURCE_FIELD
4822  equations_set_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
4823 
4824  equations_matrices=>equations%EQUATIONS_MATRICES
4825  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
4826  stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
4827  damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
4828  rhs_vector=>equations_matrices%RHS_VECTOR
4829  source_vector=>equations_matrices%SOURCE_VECTOR
4830  stiffness_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
4831  damping_matrix%ELEMENT_MATRIX%MATRIX=0.0_dp
4832  equations_mapping=>equations%EQUATIONS_MAPPING
4833 
4834 
4835  CALL field_parameter_set_data_get(equations_set_field,field_u_variable_type, &
4836  & field_values_set_type,equations_set_field_data,err,error,*999)
4837 
4838  my_compartment = equations_set_field_data(1)
4839  ncompartments = equations_set_field_data(2)
4840 
4841  linear_matrices=>equations_matrices%LINEAR_MATRICES
4842  linear_mapping=>equations_mapping%LINEAR_MAPPING
4843 
4844 
4845  num_var_count=0
4846  DO imatrix = 1,ncompartments
4847  IF(imatrix/=my_compartment)THEN
4848  num_var_count=num_var_count+1
4849  coupling_matrices(num_var_count)%PTR=>linear_matrices%MATRICES(num_var_count)%PTR
4850  field_variables(num_var_count)%PTR=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(num_var_count)%VARIABLE
4851  field_var_types(num_var_count)=field_variables(num_var_count)%PTR%VARIABLE_TYPE
4852  coupling_matrices(num_var_count)%PTR%ELEMENT_MATRIX%MATRIX=0.0_dp
4853  ENDIF
4854  END DO
4855 
4856 
4857  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
4858  field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
4859  field_var_type=field_variable%VARIABLE_TYPE
4860  geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
4861  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4862  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4863  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
4864  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4865  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
4866  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4867  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4868  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4869  & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4870  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4871  & materials_interp_parameters(field_v_variable_type)%PTR,err,error,*999)
4872  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
4873  & source_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
4874 
4875  !Loop over gauss points
4876  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
4877  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4878  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
4879  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
4880  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
4881  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4882  & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
4883  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4884  & materials_interp_point(field_v_variable_type)%PTR,err,error,*999)
4885  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
4886  & source_interp_point(field_u_variable_type)%PTR,err,error,*999)
4887  !Calculate RWG.
4888 !!TODO: Think about symmetric problems.
4889  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
4890  & quadrature_scheme%GAUSS_WEIGHTS(ng)
4891  !Loop over field components
4892  mhs=0
4893  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
4894  !Loop over element rows
4895  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4896  mhs=mhs+1
4897  nhs=0
4898  IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX) THEN
4899  !Loop over element columns
4900  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
4901  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4902  nhs=nhs+1
4903  IF(stiffness_matrix%UPDATE_MATRIX) THEN
4904  sum=0.0_dp
4905  DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
4906  pgmj(nj)=0.0_dp
4907  pgnj(nj)=0.0_dp
4908  DO ni=1,dependent_basis%NUMBER_OF_XI
4909  pgmj(nj)=pgmj(nj)+ &
4910  & quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)* &
4911  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4912  pgnj(nj)=pgnj(nj)+ &
4913  & quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)* &
4914  & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
4915  ENDDO !ni
4916  k_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR% &
4917  & values(nj,no_part_deriv)
4918  sum=sum+k_param*pgmj(nj)*pgnj(nj)
4919  ENDDO !nj
4920  coupling_param=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_v_variable_type)%PTR% &
4921  & values(my_compartment,no_part_deriv)
4922  stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+ &
4923  & sum*rwg + quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
4924  & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*rwg*coupling_param
4925  ENDIF
4926  IF(damping_matrix%UPDATE_MATRIX) THEN
4927  damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+ &
4928  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
4929  & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*rwg
4930  ENDIF
4931  ENDDO !ns
4932  ENDDO !nh
4933  ENDIF
4934  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
4935  ENDDO !ms
4936  ENDDO !mh
4937  IF(source_vector%UPDATE_VECTOR) THEN
4938  c_param=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1, no_part_deriv)
4939  mhs=0
4940  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
4941  !Loop over element rows
4942  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
4943  mhs=mhs+1
4944  source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+ &
4945  & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*c_param*rwg
4946  ENDDO !ms
4947  ENDDO !mh
4948  ENDIF
4949  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
4950  !Calculate the coupling matrices
4951 
4952  !Loop over element rows
4953  mhs=0
4954  DO mh=1,field_variable%NUMBER_OF_COMPONENTS !field_variable is the variable associated with the equations set under consideration
4955 
4956  mesh_component_1 = field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
4957  dependent_basis_1 => dependent_field%DECOMPOSITION%DOMAIN(mesh_component_1)%PTR% &
4958  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
4959  quadrature_scheme_1 => dependent_basis_1%QUADRATURE% &
4960  & quadrature_scheme_map(