OpenCMISS-Iron Internal API Documentation
Hamilton_Jacobi_equations_routines.f90
Go to the documentation of this file.
1 
43 
45 MODULE hamilton_jacobi_equations_routines
46 
47  USE base_routines
48  USE basis_routines
49  USE boundary_conditions_routines
50  USE constants
51  USE control_loop_routines
52  USE distributed_matrix_vector
53  USE domain_mappings
54  USE equations_routines
55  USE equations_mapping_routines
56  USE equations_matrices_routines
57  USE equations_set_constants
58  USE field_routines
59  USE input_output
60  USE iso_varying_string
61  USE kinds
62  USE matrix_vector
63  USE maths
64  USE node_routines
65  USE problem_constants
66  USE strings
67  USE solver_routines
68  USE timer
69  USE types
70 
71 #include "macros.h"
72 
73  IMPLICIT NONE
74 
75  PRIVATE
76 
77  !Module parameters
78 
79  !Module types
80 
81  !Module variables
82 
83  !Interfaces
84 
85 !!MERGE: move
86 
87  PUBLIC hj_boundaryconditionsanalyticcalculate
88 
89  PUBLIC hjequation_equationssetsolutionmethodset
90 
91  PUBLIC hj_equation_equations_set_setup
92 
93  PUBLIC hjequation_equationssetspecificationset
94 
95  PUBLIC hj_equation_finite_element_calculate
96 
97  PUBLIC hj_equation_problem_setup
98 
99  PUBLIC hjequation_problemspecificationset
100 
101  PUBLIC number_of_input_nodes,pre_process_information,solve_problem_fmm,solve_problem_geodesic
102  PUBLIC solve_problem_geodesic_connectivity,solve_problem_fmm_connectivity
103  PUBLIC find_minimax,post_process_data
104 
105 CONTAINS
106 
107  !
108  !================================================================================================================================
109  !
110 
111 
113  SUBROUTINE hj_boundaryconditionsanalyticcalculate(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*)
114 
115  !Argument variables
116  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
117  TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
118  INTEGER(INTG), INTENT(OUT) :: ERR
119  TYPE(varying_string), INTENT(OUT) :: ERROR
120  !Local Variables
121  INTEGER(INTG) :: component_idx,deriv_idx,dim_idx,local_ny,node_idx,NUMBER_OF_DIMENSIONS,variable_idx,variable_type
122  REAL(DP) :: VALUE,X(3)
123  REAL(DP), POINTER :: GEOMETRIC_PARAMETERS(:)
124  TYPE(domain_type), POINTER :: DOMAIN
125  TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
126  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
127  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
128  TYPE(varying_string) :: LOCAL_ERROR
129 
130  enters("HJ_BoundaryConditionsAnalyticCalculate",err,error,*999)
131 
132  IF(ASSOCIATED(equations_set)) THEN
133  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
134  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
135  IF(ASSOCIATED(dependent_field)) THEN
136  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
137  IF(ASSOCIATED(geometric_field)) THEN
138  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
139  NULLIFY(geometric_variable)
140  CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
141  CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,geometric_parameters, &
142  & err,error,*999)
143  IF(ASSOCIATED(boundary_conditions)) THEN
144  DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES
145  variable_type=dependent_field%VARIABLES(variable_idx)%VARIABLE_TYPE
146  field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
147  IF(ASSOCIATED(field_variable)) THEN
148  CALL field_parameter_set_create(dependent_field,variable_type,field_analytic_values_set_type,err,error,*999)
149  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
150  IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
151  domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
152  IF(ASSOCIATED(domain)) THEN
153  IF(ASSOCIATED(domain%TOPOLOGY)) THEN
154  domain_nodes=>domain%TOPOLOGY%NODES
155  IF(ASSOCIATED(domain_nodes)) THEN
156  !Loop over the local nodes excluding the ghosts.
157  DO node_idx=1,domain_nodes%NUMBER_OF_NODES
158  !!TODO \todo We should interpolate the geometric field here and the node position.
159  DO dim_idx=1,number_of_dimensions
160  !Default to version 1 of each node derivative
161  local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
162  & nodes(node_idx)%DERIVATIVES(1)%VERSIONS(1)
163  x(dim_idx)=geometric_parameters(local_ny)
164  ENDDO !dim_idx
165  !Loop over the derivatives
166  DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
167  SELECT CASE(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE)
168  CASE(equations_set_hj_equation_two_dim_1)
169  !u=x^2+2.x.y-y^2
170  SELECT CASE(variable_type)
171  CASE(field_u_variable_type)
172  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
173  CASE(no_global_deriv)
174  VALUE=x(1)*x(1)-2.0_dp*x(1)*x(2)-x(2)*x(2)
175  CASE(global_deriv_s1)
176  VALUE=2.0_dp*x(1)+2.0_dp*x(2)
177  CASE(global_deriv_s2)
178  VALUE=2.0_dp*x(1)-2.0_dp*x(2)
179  CASE(global_deriv_s1_s2)
180  VALUE=2.0_dp
181  CASE DEFAULT
182  local_error="The global derivative index of "//trim(number_to_vstring( &
183  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
184  & err,error))//" is invalid."
185  CALL flagerror(local_error,err,error,*999)
186  END SELECT
187  CASE(field_deludeln_variable_type)
188  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
189  CASE(no_global_deriv)
190  VALUE=0.0_dp !!TODO
191  CASE(global_deriv_s1)
192  CALL flagerror("Not implemented.",err,error,*999)
193  CASE(global_deriv_s2)
194  CALL flagerror("Not implemented.",err,error,*999)
195  CASE(global_deriv_s1_s2)
196  CALL flagerror("Not implemented.",err,error,*999)
197  CASE DEFAULT
198  local_error="The global derivative index of "//trim(number_to_vstring( &
199  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
200  & err,error))//" is invalid."
201  CALL flagerror(local_error,err,error,*999)
202  END SELECT
203  CASE DEFAULT
204  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
205  & " is invalid."
206  CALL flagerror(local_error,err,error,*999)
207  END SELECT
208  CASE(equations_set_hj_equation_two_dim_2)
209  !u=cos(x).cosh(y)
210  SELECT CASE(variable_type)
211  CASE(field_u_variable_type)
212  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
213  CASE(no_global_deriv)
214  VALUE=cos(x(1))*cosh(x(2))
215  CASE(global_deriv_s1)
216  VALUE=-sin(x(1))*cosh(x(2))
217  CASE(global_deriv_s2)
218  VALUE=cos(x(1))*sinh(x(2))
219  CASE(global_deriv_s1_s2)
220  VALUE=-sin(x(1))*sinh(x(2))
221  CASE DEFAULT
222  local_error="The global derivative index of "//trim(number_to_vstring( &
223  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
224  & err,error))//" is invalid."
225  CALL flagerror(local_error,err,error,*999)
226  END SELECT
227  CASE(field_deludeln_variable_type)
228  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
229  CASE(no_global_deriv)
230  VALUE=0.0_dp !!TODO
231  CASE(global_deriv_s1)
232  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
233  CASE(global_deriv_s2)
234  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
235  CASE(global_deriv_s1_s2)
236  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
237  CASE DEFAULT
238  local_error="The global derivative index of "//trim(number_to_vstring( &
239  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
240  & err,error))//" is invalid."
241  CALL flagerror(local_error,err,error,*999)
242  END SELECT
243  CASE DEFAULT
244  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
245  & " is invalid."
246  CALL flagerror(local_error,err,error,*999)
247  END SELECT
248  CASE(equations_set_hj_equation_three_dim_1)
249  !u=x^2+y^2-2.z^2
250  SELECT CASE(variable_type)
251  CASE(field_u_variable_type)
252  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
253  CASE(no_global_deriv)
254  VALUE=x(1)*x(1)+x(2)*x(2)-2.0_dp*x(3)*x(3)
255  CASE(global_deriv_s1)
256  VALUE=2.0_dp*x(1)
257  CASE(global_deriv_s2)
258  VALUE=2.0_dp*x(2)
259  CASE(global_deriv_s1_s2)
260  VALUE=0.0_dp
261  CASE(global_deriv_s3)
262  VALUE=-4.0_dp*x(3)
263  CASE(global_deriv_s1_s3)
264  VALUE=0.0_dp
265  CASE(global_deriv_s2_s3)
266  VALUE=0.0_dp
267  CASE(global_deriv_s1_s2_s3)
268  VALUE=0.0_dp
269  CASE DEFAULT
270  local_error="The global derivative index of "//trim(number_to_vstring( &
271  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
272  & err,error))//" is invalid."
273  CALL flagerror(local_error,err,error,*999)
274  END SELECT
275  CASE(field_deludeln_variable_type)
276  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
277  CASE(no_global_deriv)
278  VALUE=0.0_dp !!TODO
279  CASE(global_deriv_s1)
280  CALL flagerror("Not implemented.",err,error,*999)
281  CASE(global_deriv_s2)
282  CALL flagerror("Not implemented.",err,error,*999)
283  CASE(global_deriv_s1_s2)
284  CALL flagerror("Not implemented.",err,error,*999)
285  CASE(global_deriv_s3)
286  CALL flagerror("Not implemented.",err,error,*999)
287  CASE(global_deriv_s1_s3)
288  CALL flagerror("Not implemented.",err,error,*999)
289  CASE(global_deriv_s2_s3)
290  CALL flagerror("Not implemented.",err,error,*999)
291  CASE(global_deriv_s1_s2_s3)
292  CALL flagerror("Not implemented.",err,error,*999)
293  CASE DEFAULT
294  local_error="The global derivative index of "//trim(number_to_vstring( &
295  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
296  & err,error))//" is invalid."
297  CALL flagerror(local_error,err,error,*999)
298  END SELECT
299  CASE DEFAULT
300  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
301  & " is invalid."
302  CALL flagerror(local_error,err,error,*999)
303  END SELECT
304  CASE(equations_set_hj_equation_three_dim_2)
305  !u=cos(x).cosh(y).z
306  SELECT CASE(variable_type)
307  CASE(field_u_variable_type)
308  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
309  CASE(no_global_deriv)
310  VALUE=cos(x(1))*cosh(x(2))*x(3)
311  CASE(global_deriv_s1)
312  VALUE=-sin(x(1))*cosh(x(2))*x(3)
313  CASE(global_deriv_s2)
314  VALUE=cos(x(1))*sinh(x(2))*x(3)
315  CASE(global_deriv_s1_s2)
316  VALUE=-sin(x(1))*sinh(x(2))*x(3)
317  CASE(global_deriv_s3)
318  VALUE=cos(x(1))*cosh(x(2))
319  CASE(global_deriv_s1_s3)
320  VALUE=-sin(x(1))*cosh(x(2))
321  CASE(global_deriv_s2_s3)
322  VALUE=cos(x(1))*sinh(x(2))
323  CASE(global_deriv_s1_s2_s3)
324  VALUE=-sin(x(1))*sinh(x(2))
325  CASE DEFAULT
326  local_error="The global derivative index of "//trim(number_to_vstring( &
327  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
328  & err,error))//" is invalid."
329  CALL flagerror(local_error,err,error,*999)
330  END SELECT
331  CASE(field_deludeln_variable_type)
332  SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
333  CASE(no_global_deriv)
334  VALUE=0.0_dp !!TODO
335  CASE(global_deriv_s1)
336  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
337  CASE(global_deriv_s2)
338  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
339  CASE(global_deriv_s1_s2)
340  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
341  CASE(global_deriv_s3)
342  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
343  CASE(global_deriv_s1_s3)
344  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
345  CASE(global_deriv_s2_s3)
346  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
347  CASE(global_deriv_s1_s2_s3)
348  !CALL FlagError("Not implemented.",ERR,ERROR,*999)
349  CASE DEFAULT
350  local_error="The global derivative index of "//trim(number_to_vstring( &
351  domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
352  & err,error))//" is invalid."
353  CALL flagerror(local_error,err,error,*999)
354  END SELECT
355  CASE DEFAULT
356  local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))// &
357  & " is invalid."
358  CALL flagerror(local_error,err,error,*999)
359  END SELECT
360  CASE DEFAULT
361  local_error="The analytic function type of "// &
362  & trim(number_to_vstring(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
363  & " is invalid."
364  CALL flagerror(local_error,err,error,*999)
365  END SELECT
366  !Default to version 1 of each node derivative
367  local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
368  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
369  CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
370  & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
371  IF(variable_type==field_u_variable_type) THEN
372  IF(domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN
373  !If we are a boundary node then set the analytic value on the boundary
374  CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field,variable_type, &
375  & local_ny,boundary_condition_fixed,VALUE,err,error,*999)
376  ENDIF
377  ENDIF
378  ENDDO !deriv_idx
379  ENDDO !node_idx
380  ELSE
381  CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
382  ENDIF
383  ELSE
384  CALL flagerror("Domain topology is not associated.",err,error,*999)
385  ENDIF
386  ELSE
387  CALL flagerror("Domain is not associated.",err,error,*999)
388  ENDIF
389  ELSE
390  CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
391  ENDIF
392  ENDDO !component_idx
393  CALL field_parameter_set_update_start(dependent_field,variable_type,field_analytic_values_set_type, &
394  & err,error,*999)
395  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_analytic_values_set_type, &
396  & err,error,*999)
397  ELSE
398  CALL flagerror("Field variable is not associated.",err,error,*999)
399  ENDIF
400 
401  ENDDO !variable_idx
402  CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
403  & geometric_parameters,err,error,*999)
404  ELSE
405  CALL flagerror("Boundary conditions is not associated.",err,error,*999)
406  ENDIF
407  ELSE
408  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
409  ENDIF
410  ELSE
411  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
412  ENDIF
413  ELSE
414  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
415  ENDIF
416  ELSE
417  CALL flagerror("Equations set is not associated.",err,error,*999)
418  ENDIF
419 
420  exits("HJ_BoundaryConditionsAnalyticCalculate")
421  RETURN
422 999 errorsexits("HJ_BoundaryConditionsAnalyticCalculate",err,error)
423  RETURN 1
424  END SUBROUTINE hj_boundaryconditionsanalyticcalculate
425 
426  !
427  !================================================================================================================================
428  !
429 
431  SUBROUTINE hj_equation_finite_element_calculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
432 
433  !Argument variables
434  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
435  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
436  INTEGER(INTG), INTENT(OUT) :: ERR
437  TYPE(varying_string), INTENT(OUT) :: ERROR
438  !Local Variables
439  INTEGER(INTG) FIELD_VAR_TYPE,ng,mh,mhs,mi,ms,nh,nhs,ni,ns
440  REAL(DP) :: RWG,SUM,PGMSI(3),PGNSI(3)
441  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
442  TYPE(equations_type), POINTER :: EQUATIONS
443  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
444  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
445  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
446  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
447  TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
448  TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
449  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
450  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
451  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
452  TYPE(varying_string) :: LOCAL_ERROR
453 
454 #ifdef TAUPROF
455  CHARACTER(26) :: CVAR
456  INTEGER :: GAUSS_POINT_LOOP_PHASE(2) = (/ 0, 0 /)
457  SAVE gauss_point_loop_phase
458 #endif
459 
460  enters("HJ_EQUATION_FINITE_ELEMENT_CALCULATE",err,error,*999)
461 
462  IF(ASSOCIATED(equations_set)) THEN
463  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
464  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
465  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
466  CALL flagerror("Equations set specification must have three entries for a Hamilton-Jacobi type equations set.", &
467  & err,error,*999)
468  END IF
469  equations=>equations_set%EQUATIONS
470  IF(ASSOCIATED(equations)) THEN
471  SELECT CASE(equations_set%SPECIFICATION(3))
472  CASE(equations_set_standard_hj_subtype)
473 !!TODO: move these and scale factor adjustment out once generalised Hamilton-Jacobi is put in.
474  !Store all these in equations matrices/somewhere else?????
475  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
476  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
477  equations_matrices=>equations%EQUATIONS_MATRICES
478  linear_matrices=>equations_matrices%LINEAR_MATRICES
479  equations_matrix=>linear_matrices%MATRICES(1)%PTR
480  rhs_vector=>equations_matrices%RHS_VECTOR
481  equations_mapping=>equations%EQUATIONS_MAPPING
482  linear_mapping=>equations_mapping%LINEAR_MAPPING
483  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
484  field_var_type=field_variable%VARIABLE_TYPE
485  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
486  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
487  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
488  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
489  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
490  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
491  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
492  !Loop over gauss points
493  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
494 #ifdef TAUPROF
495  WRITE (cvar,'(a17,i2)') 'Gauss Point Loop ',ng
496  CALL tau_phase_create_dynamic(gauss_point_loop_phase,cvar)
497  CALL tau_phase_start(gauss_point_loop_phase)
498 #endif
499  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
500  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
501  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
502  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
503  !Calculate RWG.
504 !!TODO: Think about symmetric problems.
505  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
506  & quadrature_scheme%GAUSS_WEIGHTS(ng)
507  !Loop over field components
508  mhs=0
509  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
510  !Loop over element rows
511 !!TODO: CHANGE ELEMENT CALCULATE TO WORK OF ns ???
512  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
513  mhs=mhs+1
514  nhs=0
515  IF(equations_matrix%UPDATE_MATRIX) THEN
516  !Loop over element columns
517  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
518  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
519  nhs=nhs+1
520  DO ni=1,dependent_basis%NUMBER_OF_XI
521  pgmsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
522  pgnsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
523  ENDDO !ni
524 
525  sum=0.0_dp
526  DO mi=1,dependent_basis%NUMBER_OF_XI
527  DO ni=1,dependent_basis%NUMBER_OF_XI
528  sum=sum+pgmsi(mi)*pgnsi(ni)*equations%INTERPOLATION% &
529  & geometric_interp_point_metrics(field_u_variable_type)%PTR%GU(mi,ni)
530  ENDDO !ni
531  ENDDO !mi
532  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg
533 
534  ENDDO !ns
535  ENDDO !nh
536  ENDIF
537  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
538  ENDDO !ms
539  ENDDO !mh
540 #ifdef TAUPROF
541  CALL tau_phase_stop(gauss_point_loop_phase)
542 #endif
543  ENDDO !ng
544 
545  !Scale factor adjustment
546  IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
547  CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
548  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
549  mhs=0
550  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
551  !Loop over element rows
552  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
553  mhs=mhs+1
554  nhs=0
555  IF(equations_matrix%UPDATE_MATRIX) THEN
556  !Loop over element columns
557  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
558  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
559  nhs=nhs+1
560  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
561  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
562  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
563  ENDDO !ns
564  ENDDO !nh
565  ENDIF
566  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
567  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
568  ENDDO !ms
569  ENDDO !mh
570  ENDIF
571  CASE(equations_set_generalised_hj_subtype)
572  CALL flagerror("Not implemented.",err,error,*999)
573  CASE DEFAULT
574  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
575  & " is not valid for a Hamilton-Jacobi equation type of a classical field equations set class."
576  CALL flagerror(local_error,err,error,*999)
577  END SELECT
578 
579  ELSE
580  CALL flagerror("Equations set equations is not associated.",err,error,*999)
581  ENDIF
582  ELSE
583  CALL flagerror("Equations set is not associated.",err,error,*999)
584  ENDIF
585 
586  exits("HJ_EQUATION_FINITE_ELEMENT_CALCULATE")
587  RETURN
588 999 errorsexits("HJ_EQUATION_FINITE_ELEMENT_CALCULATE",err,error)
589  RETURN 1
590  END SUBROUTINE hj_equation_finite_element_calculate
591 
592  !
593  !================================================================================================================================
594  !
595 
597  SUBROUTINE hj_equation_fast_marching_calculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
598 
599  !Argument variables
600  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
601  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
602  INTEGER(INTG), INTENT(OUT) :: ERR
603  TYPE(varying_string), INTENT(OUT) :: ERROR
604  !Local Variables
605  INTEGER(INTG) FIELD_VAR_TYPE,ng,mh,mhs,mi,ms,nh,nhs,ni,ns
606  REAL(DP) :: RWG,SUM,PGMSI(3),PGNSI(3)
607  TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
608  TYPE(equations_type), POINTER :: EQUATIONS
609  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
610  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
611  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
612  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
613  TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
614  TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
615  TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
616  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
617  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
618  TYPE(varying_string) :: LOCAL_ERROR
619 
620 #ifdef TAUPROF
621  CHARACTER(26) :: CVAR
622  INTEGER :: GAUSS_POINT_LOOP_PHASE(2) = (/ 0, 0 /)
623  SAVE gauss_point_loop_phase
624 #endif
625 
626  enters("HJ_EQUATION_FAST_MARCHING_CALCULATE",err,error,*999)
627 
628  IF(ASSOCIATED(equations_set)) THEN
629  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
630  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
631  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
632  CALL flagerror("Equations set specification must have three entries for a Hamilton-Jacobi type equations set.", &
633  & err,error,*999)
634  END IF
635  equations=>equations_set%EQUATIONS
636  IF(ASSOCIATED(equations)) THEN
637  SELECT CASE(equations_set%SPECIFICATION(3))
638  CASE(equations_set_standard_hj_subtype)
639 !!TODO: move these and scale factor adjustment out once generalised Hamilton-Jacobi is put in.
640  !Store all these in equations matrices/somewhere else?????
641  dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
642 ! MATERIALS_FIELD=>EQUATIONS%INTERPOLATION%MATERIALS_FIELD!
643  geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
644  equations_matrices=>equations%EQUATIONS_MATRICES
645  linear_matrices=>equations_matrices%LINEAR_MATRICES
646  equations_matrix=>linear_matrices%MATRICES(1)%PTR
647  rhs_vector=>equations_matrices%RHS_VECTOR
648  equations_mapping=>equations%EQUATIONS_MAPPING
649  linear_mapping=>equations_mapping%LINEAR_MAPPING
650  field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
651  field_var_type=field_variable%VARIABLE_TYPE
652  dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
653  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
654  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
655  & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
656  quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
657  CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
658  & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
659  !Loop over gauss points
660  DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
661 #ifdef TAUPROF
662  WRITE (cvar,'(a17,i2)') 'Gauss Point Loop ',ng
663  CALL tau_phase_create_dynamic(gauss_point_loop_phase,cvar)
664  CALL tau_phase_start(gauss_point_loop_phase)
665 #endif
666  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
667  & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
668  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
669  & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
670  !Calculate RWG.
671 !!TODO: Think about symmetric problems.
672  rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
673  & quadrature_scheme%GAUSS_WEIGHTS(ng)
674  !Loop over field components
675  mhs=0
676  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
677  !Loop over element rows
678 !!TODO: CHANGE ELEMENT CALCULATE TO WORK OF ns ???
679  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
680  mhs=mhs+1
681  nhs=0
682  IF(equations_matrix%UPDATE_MATRIX) THEN
683  !Loop over element columns
684  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
685  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
686  nhs=nhs+1
687  DO ni=1,dependent_basis%NUMBER_OF_XI
688  pgmsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
689  pgnsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
690  ENDDO !ni
691 
692  sum=0.0_dp
693  DO mi=1,dependent_basis%NUMBER_OF_XI
694  DO ni=1,dependent_basis%NUMBER_OF_XI
695  sum=sum+pgmsi(mi)*pgnsi(ni)*equations%INTERPOLATION% &
696  & geometric_interp_point_metrics(field_u_variable_type)%PTR%GU(mi,ni)
697  ENDDO !ni
698  ENDDO !mi
699  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg
700 
701  ENDDO !ns
702  ENDDO !nh
703  ENDIF
704  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
705  ENDDO !ms
706  ENDDO !mh
707 #ifdef TAUPROF
708  CALL tau_phase_stop(gauss_point_loop_phase)
709 #endif
710  ENDDO !ng
711 
712  !Scale factor adjustment
713  IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
714  CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
715  & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
716  mhs=0
717  DO mh=1,field_variable%NUMBER_OF_COMPONENTS
718  !Loop over element rows
719  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
720  mhs=mhs+1
721  nhs=0
722  IF(equations_matrix%UPDATE_MATRIX) THEN
723  !Loop over element columns
724  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
725  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
726  nhs=nhs+1
727  equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
728  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
729  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
730  ENDDO !ns
731  ENDDO !nh
732  ENDIF
733  IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
734  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
735  ENDDO !ms
736  ENDDO !mh
737  ENDIF
738  CASE(equations_set_generalised_hj_subtype)
739  CALL flagerror("Not implemented.",err,error,*999)
740  CASE DEFAULT
741  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
742  & " is not valid for a Hamilton-Jacobi equation type of a classical field equations set class."
743  CALL flagerror(local_error,err,error,*999)
744  END SELECT
745 
746  ELSE
747  CALL flagerror("Equations set equations is not associated.",err,error,*999)
748  ENDIF
749  ELSE
750  CALL flagerror("Equations set is not associated.",err,error,*999)
751  ENDIF
752 
753  exits("HJ_EQUATION_FAST_MARCHING_CALCULATE")
754  RETURN
755 999 errorsexits("HJ_EQUATION_FAST_MARCHING_CALCULATE",err,error)
756  RETURN 1
757  END SUBROUTINE hj_equation_fast_marching_calculate
758 
759  !
760  !================================================================================================================================
761  !
762 
764  SUBROUTINE hj_equation_equations_set_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
765 
766  !Argument variables
767  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
768  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
769  INTEGER(INTG), INTENT(OUT) :: ERR
770  TYPE(varying_string), INTENT(OUT) :: ERROR
771  !Local Variables
772  TYPE(varying_string) :: LOCAL_ERROR
773 
774  enters("HJ_EQUATION_EQUATIONS_SET_SETUP",err,error,*999)
775 
776  IF(ASSOCIATED(equations_set)) THEN
777  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
778  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
779  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
780  CALL flagerror("Equations set specification must have three entries for a Hamilton-Jacobi type equations set.", &
781  & err,error,*999)
782  END IF
783  SELECT CASE(equations_set%SPECIFICATION(3))
784 
785  CASE(equations_set_standard_hj_subtype)
786  CALL hj_equation_equations_set_standard_setup(equations_set,equations_set_setup,err,error,*999)
787  CASE(equations_set_generalised_hj_subtype)
788  CALL flagerror("Not implemented.",err,error,*999)
789  CASE DEFAULT
790  local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
791  & " is not valid for a Hamilton-Jacobi equation type of a classical field equation set class."
792  CALL flagerror(local_error,err,error,*999)
793  END SELECT
794  ELSE
795  CALL flagerror("Equations set is not associated.",err,error,*999)
796  ENDIF
797 
798  exits("HJ_EQUATION_EQUATIONS_SET_SETUP")
799  RETURN
800 999 errorsexits("HJ_EQUATION_EQUATIONS_SET_SETUP",err,error)
801  RETURN 1
802  END SUBROUTINE hj_equation_equations_set_setup
803 
804  !
805  !================================================================================================================================
806  !
807 
809  SUBROUTINE hjequation_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)
810 
811  !Argument variables
812  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
813  INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
814  INTEGER(INTG), INTENT(OUT) :: ERR
815  TYPE(varying_string), INTENT(OUT) :: ERROR
816  !Local Variables
817  TYPE(varying_string) :: LOCAL_ERROR
818 
819  enters("HJEquation_EquationsSetSolutionMethodSet",err,error,*999)
820 
821  IF(ASSOCIATED(equations_set)) THEN
822  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
823  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
824  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
825  CALL flagerror("Equations set specification must have three entries for a Hamilton-Jacobi type equations set.", &
826  & err,error,*999)
827  END IF
828  SELECT CASE(equations_set%SPECIFICATION(3))
829  CASE(equations_set_standard_hj_subtype)
830  SELECT CASE(solution_method)
831  CASE(equations_set_fem_solution_method)
832  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
833 ! CASE(EQUATIONS_SET_FMM_SOLUTION_METHOD)
834 ! CALL FlagError("Not implemented.",ERR,ERROR,*999)
835  CASE(equations_set_bem_solution_method)
836  CALL flagerror("Not implemented.",err,error,*999)
837  CASE(equations_set_fd_solution_method)
838  CALL flagerror("Not implemented.",err,error,*999)
839  CASE(equations_set_fv_solution_method)
840  CALL flagerror("Not implemented.",err,error,*999)
841  CASE(equations_set_gfem_solution_method)
842  CALL flagerror("Not implemented.",err,error,*999)
843  CASE(equations_set_gfv_solution_method)
844  CALL flagerror("Not implemented.",err,error,*999)
845  CASE DEFAULT
846  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
847  CALL flagerror(local_error,err,error,*999)
848  END SELECT
849  CASE(equations_set_generalised_hj_subtype)
850  SELECT CASE(solution_method)
851  CASE(equations_set_fem_solution_method)
852  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
853 ! CASE(EQUATIONS_SET_FMM_SOLUTION_METHOD)
854 ! CALL FlagError("Not implemented.",ERR,ERROR,*999)
855  CASE(equations_set_bem_solution_method)
856  CALL flagerror("Not implemented.",err,error,*999)
857  CASE(equations_set_fd_solution_method)
858  CALL flagerror("Not implemented.",err,error,*999)
859  CASE(equations_set_fv_solution_method)
860  CALL flagerror("Not implemented.",err,error,*999)
861  CASE(equations_set_gfem_solution_method)
862  CALL flagerror("Not implemented.",err,error,*999)
863  CASE(equations_set_gfv_solution_method)
864  CALL flagerror("Not implemented.",err,error,*999)
865  CASE DEFAULT
866  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
867  CALL flagerror(local_error,err,error,*999)
868  END SELECT
869  CASE DEFAULT
870  local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
871  & " is not valid for a Hamilton-Jacobi equation type of an classical field equations set class."
872  CALL flagerror(local_error,err,error,*999)
873  END SELECT
874  ELSE
875  CALL flagerror("Equations set is not associated.",err,error,*999)
876  ENDIF
877 
878  exits("HJEquation_EquationsSetSolutionMethodSet")
879  RETURN
880 999 errorsexits("HJEquation_EquationsSetSolutionMethodSet",err,error)
881  RETURN 1
882  END SUBROUTINE hjequation_equationssetsolutionmethodset
883 
884  !
885  !================================================================================================================================
886  !
887 
889  SUBROUTINE hjequation_equationssetspecificationset(equationsSet,specification,err,error,*)
890 
891  !Argument variables
892  TYPE(equations_set_type), POINTER :: equationsSet
893  INTEGER(INTG), INTENT(IN) :: specification(:)
894  INTEGER(INTG), INTENT(OUT) :: err
895  TYPE(varying_string), INTENT(OUT) :: error
896  !Local Variables
897  TYPE(varying_string) :: localError
898  INTEGER(INTG) :: subtype
899 
900  CALL enters("HJEquation_EquationsSetSpecificationSet",err,error,*999)
901 
902  IF(ASSOCIATED(equationsset)) THEN
903  IF(SIZE(specification,1)/=3) THEN
904  CALL flagerror("Equations set specification must have three entries for a Hamilton-Jacobi type equations set.", &
905  & err,error,*999)
906  END IF
907  subtype=specification(3)
908  SELECT CASE(subtype)
909  CASE(equations_set_standard_hj_subtype)
910  !ok
911  CASE(equations_set_generalised_hj_subtype)
912  CALL flagerror("Not implemented.",err,error,*999)
913  CASE DEFAULT
914  localerror="The third equations set specification of "//trim(numbertovstring(subtype,"*",err,error))// &
915  & " is not valid for a Hamilton-Jacobi type of a classical field equations set."
916  CALL flagerror(localerror,err,error,*999)
917  END SELECT
918  !Set full specification
919  IF(ALLOCATED(equationsset%specification)) THEN
920  CALL flagerror("Equations set specification is already allocated.",err,error,*999)
921  ELSE
922  ALLOCATE(equationsset%specification(3),stat=err)
923  IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
924  END IF
925  equationsset%specification(1:3)=[equations_set_classical_field_class,equations_set_hj_equation_type,subtype]
926  ELSE
927  CALL flagerror("Equations set is not associated.",err,error,*999)
928  END IF
929 
930  exits("HJEquation_EquationsSetSpecificationSet")
931  RETURN
932 999 errors("HJEquation_EquationsSetSpecificationSet",err,error)
933  exits("HJEquation_EquationsSetSpecificationSet")
934  RETURN 1
935 
936  END SUBROUTINE hjequation_equationssetspecificationset
937 
938  !
939  !================================================================================================================================
940  !
941 
943  SUBROUTINE hj_equation_equations_set_standard_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
944 
945  !Argument variables
946  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
947  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
948  INTEGER(INTG), INTENT(OUT) :: ERR
949  TYPE(varying_string), INTENT(OUT) :: ERROR
950  !Local Variables
951  INTEGER(INTG) :: component_idx,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE,NUMBER_OF_DIMENSIONS, &
952  & NUMBER_OF_MATERIALS_COMPONENTS, GEOMETRIC_MESH_COMPONENT
953  TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
954  TYPE(field_type), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,GEOMETRIC_FIELD
955  TYPE(equations_type), POINTER :: EQUATIONS
956  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
957  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
958  TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
959 
960  TYPE(varying_string) :: LOCAL_ERROR
961 
962  enters("HJ_EQUATION_EQUATION_SET_STANDARD_SETUP",err,error,*999)
963 
964  NULLIFY(equations)
965  NULLIFY(equations_mapping)
966  NULLIFY(equations_matrices)
967  NULLIFY(geometric_decomposition)
968 
969  IF(ASSOCIATED(equations_set)) THEN
970  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
971  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
972  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
973  CALL flagerror("Equations set specification must have three entries for a Hamilton-Jacobi type equations set.", &
974  & err,error,*999)
975  END IF
976  IF(equations_set%SPECIFICATION(3)==equations_set_standard_hj_subtype) THEN
977  SELECT CASE(equations_set_setup%SETUP_TYPE)
978  CASE(equations_set_setup_initial_type)
979  SELECT CASE(equations_set_setup%ACTION_TYPE)
980  CASE(equations_set_setup_start_action)
981  CALL hjequation_equationssetsolutionmethodset(equations_set,equations_set_fem_solution_method,err,error,*999)
982  CASE(equations_set_setup_finish_action)
983  !Do nothing
984  CASE DEFAULT
985  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
986  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
987  & " is invalid for a standard Hamilton-Jacobi equation."
988  CALL flagerror(local_error,err,error,*999)
989  END SELECT
990  CASE(equations_set_setup_geometry_type)
991  !Do nothing
992  CASE(equations_set_setup_dependent_type)
993  SELECT CASE(equations_set_setup%ACTION_TYPE)
994  CASE(equations_set_setup_start_action)
995  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
996  !Create the auto created dependent field
997  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
998  & dependent_field,err,error,*999)
999  CALL field_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,"Dependent Field",err,error,*999)
1000  CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
1001  CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
1002  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1003  CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
1004  & err,error,*999)
1005  CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
1006  & geometric_field,err,error,*999)
1007  CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
1008  CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,(/field_u_variable_type, &
1009  & field_deludeln_variable_type/),err,error,*999)
1010  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,"Phi",err,error,*999)
1011  CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,"del Phi/del n", &
1012  & err,error,*999)
1013  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
1014  & field_scalar_dimension_type,err,error,*999)
1015  CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
1016  & field_scalar_dimension_type,err,error,*999)
1017  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
1018  & field_dp_type,err,error,*999)
1019  CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
1020  & field_dp_type,err,error,*999)
1021  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
1022  & err,error,*999)
1023  CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
1024  & err,error,*999)
1025  CALL field_component_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1,"Phi",err,error,*999)
1026  CALL field_component_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
1027  & "del Phi/del n",err,error,*999)
1028  !Default to the geometric interpolation setup
1029  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
1030  & geometric_mesh_component,err,error,*999)
1031  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
1032  & geometric_mesh_component,err,error,*999)
1033  CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
1034  & geometric_mesh_component,err,error,*999)
1035  SELECT CASE(equations_set%SOLUTION_METHOD)
1036  CASE(equations_set_fem_solution_method)
1037  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
1038  & field_u_variable_type,1,field_node_based_interpolation,err,error,*999)
1039  CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
1040  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
1041  !Default the scaling to the geometric field scaling
1042  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1043  CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
1044 ! CASE(EQUATIONS_SET_FMM_SOLUTION_METHOD)
1045 ! CALL FlagError("Not implemented.",ERR,ERROR,*999)
1046  CASE(equations_set_bem_solution_method)
1047  CALL flagerror("Not implemented.",err,error,*999)
1048  CASE(equations_set_fd_solution_method)
1049  CALL flagerror("Not implemented.",err,error,*999)
1050  CASE(equations_set_fv_solution_method)
1051  CALL flagerror("Not implemented.",err,error,*999)
1052  CASE(equations_set_gfem_solution_method)
1053  CALL flagerror("Not implemented.",err,error,*999)
1054  CASE(equations_set_gfv_solution_method)
1055  CALL flagerror("Not implemented.",err,error,*999)
1056  CASE DEFAULT
1057  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1058  & " is invalid."
1059  CALL flagerror(local_error,err,error,*999)
1060  END SELECT
1061  ELSE
1062  !Check the user specified field
1063  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
1064  CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
1065  CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
1066  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type,field_deludeln_variable_type/), &
1067  & err,error,*999)
1068  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type,err,error,*999)
1069  CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
1070  & err,error,*999)
1071  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1072  CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
1073  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
1074  CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1,err,error,*999)
1075  SELECT CASE(equations_set%SOLUTION_METHOD)
1076  CASE(equations_set_fem_solution_method)
1077  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
1078  & field_node_based_interpolation,err,error,*999)
1079  CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
1080  & field_node_based_interpolation,err,error,*999)
1081 ! CASE(EQUATIONS_SET_FMM_SOLUTION_METHOD)
1082 ! CALL FlagError("Not implemented.",ERR,ERROR,*999)
1083  CASE(equations_set_bem_solution_method)
1084  CALL flagerror("Not implemented.",err,error,*999)
1085  CASE(equations_set_fd_solution_method)
1086  CALL flagerror("Not implemented.",err,error,*999)
1087  CASE(equations_set_fv_solution_method)
1088  CALL flagerror("Not implemented.",err,error,*999)
1089  CASE(equations_set_gfem_solution_method)
1090  CALL flagerror("Not implemented.",err,error,*999)
1091  CASE(equations_set_gfv_solution_method)
1092  CALL flagerror("Not implemented.",err,error,*999)
1093  CASE DEFAULT
1094  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1095  & " is invalid."
1096  CALL flagerror(local_error,err,error,*999)
1097  END SELECT
1098  ENDIF
1099  CASE(equations_set_setup_finish_action)
1100  IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
1101  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
1102  ENDIF
1103  CASE DEFAULT
1104  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1105  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1106  & " is invalid for a standard Hamilton-Jacobi equation"
1107  CALL flagerror(local_error,err,error,*999)
1108  END SELECT
1109  CASE(equations_set_setup_materials_type)
1110  SELECT CASE(equations_set_setup%ACTION_TYPE)
1111  CASE(equations_set_setup_start_action)
1112  !Do nothing
1113 
1114 
1115 
1116 
1117  equations_materials=>equations_set%MATERIALS
1118  IF(ASSOCIATED(equations_materials)) THEN
1119  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
1120  !Create the auto created materials field
1121  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
1122  & materials_field,err,error,*999)
1123  CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
1124  CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
1125  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
1126  CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
1127  & err,error,*999)
1128  CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
1129  & geometric_field,err,error,*999)
1130  CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
1131  CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,(/field_u_variable_type/), &
1132  & err,error,*999)
1133  CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1134  & field_vector_dimension_type,err,error,*999)
1135  CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1136  & field_dp_type,err,error,*999)
1137  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1138  & number_of_dimensions,err,error,*999)
1139  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
1140  !Constant source. Materials field components are 1 for each dimension and 1 for the constant source
1141  !i.e., k and c in div(k.grad(u(x)))=c(x)
1142  number_of_materials_components=number_of_dimensions+1
1143  ELSE
1144  !Linear source. Materials field components are 1 for each dimension and 2 for the linear source
1145  !i.e., k and a and c in div(k.grad(u(x)))=a(x)u(x)+c(x)
1146  number_of_materials_components=number_of_dimensions+2
1147  ENDIF
1148  !Set the number of materials components
1149  CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1150  & number_of_materials_components,err,error,*999)
1151  !Default the k materials components to the geometric interpolation setup with constant interpolation
1152  DO component_idx=1,number_of_dimensions
1153  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1154  & component_idx,geometric_component_number,err,error,*999)
1155  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1156  & component_idx,geometric_component_number,err,error,*999)
1157  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1158  & component_idx,field_constant_interpolation,err,error,*999)
1159  ENDDO !component_idx
1160  !Default the source materials components to the first component geometric interpolation with constant interpolation
1161  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1162  & 1,geometric_component_number,err,error,*999)
1163  DO component_idx=number_of_dimensions+1,number_of_materials_components
1164  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1165  & component_idx,geometric_component_number,err,error,*999)
1166  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1167  & component_idx,field_constant_interpolation,err,error,*999)
1168  ENDDO !components_idx
1169  !Default the field scaling to that of the geometric field
1170  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
1171  CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
1172  ELSE
1173  !Check the user specified field
1174  CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
1175  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
1176  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
1177  CALL field_variable_types_check(equations_set_setup%FIELD,(/field_u_variable_type/),err,error,*999)
1178  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
1179  & err,error,*999)
1180  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
1181  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1182  & number_of_dimensions,err,error,*999)
1183  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
1184  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+1, &
1185  & err,error,*999)
1186  ELSE
1187  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+2, &
1188  & err,error,*999)
1189  ENDIF
1190  ENDIF
1191  ELSE
1192  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1193  ENDIF
1194 
1195 
1196 
1197 
1198  CASE(equations_set_setup_finish_action)
1199  !Do nothing
1200 
1201 
1202  equations_materials=>equations_set%MATERIALS
1203  IF(ASSOCIATED(equations_materials)) THEN
1204  IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
1205  !Finish creating the materials field
1206  CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
1207  !Set the default values for the materials field
1208  CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
1209  & number_of_dimensions,err,error,*999)
1210  IF(equations_set%SPECIFICATION(3)==equations_set_constant_source_poisson_subtype) THEN
1211  !Constant source. Materials field components are 1 for each dimension and 1 for the constant source
1212  !i.e., k and c in div(k.grad(u(x)))=c(x)
1213  number_of_materials_components=number_of_dimensions+1
1214  ELSE
1215  !Linear source. Materials field components are 1 for each dimension and 2 for the linear source
1216  !i.e., k and a and c in div(k.grad(u(x)))=a(x)u(x)+c(x)
1217  number_of_materials_components=number_of_dimensions+2
1218  ENDIF
1219  !First set the k values to 1.0
1220  DO component_idx=1,number_of_dimensions
1221  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1222  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
1223  ENDDO !component_idx
1224  !Now set the source values to 1.0
1225  DO component_idx=number_of_dimensions+1,number_of_materials_components
1226  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
1227  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
1228  ENDDO !component_idx
1229  ENDIF
1230  ELSE
1231  CALL flagerror("Equations set materials is not associated.",err,error,*999)
1232  ENDIF
1233 
1234 
1235 
1236  CASE DEFAULT
1237  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1238  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1239  & " is invalid for a standard Hamilton-Jacobi equation."
1240  CALL flagerror(local_error,err,error,*999)
1241  END SELECT
1242  CASE(equations_set_setup_source_type)
1243  SELECT CASE(equations_set_setup%ACTION_TYPE)
1244  CASE(equations_set_setup_start_action)
1245  !Do nothing
1246  CASE(equations_set_setup_finish_action)
1247  !Do nothing
1248  CASE DEFAULT
1249  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1250  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1251  & " is invalid for a standard Hamilton-Jacobi equation."
1252  CALL flagerror(local_error,err,error,*999)
1253  END SELECT
1254  CASE(equations_set_setup_analytic_type)
1255  SELECT CASE(equations_set_setup%ACTION_TYPE)
1256  CASE(equations_set_setup_start_action)
1257  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
1258  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
1259  IF(ASSOCIATED(dependent_field)) THEN
1260  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
1261  IF(ASSOCIATED(geometric_field)) THEN
1262  CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
1263  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
1264  CASE(equations_set_hj_equation_two_dim_1)
1265  !Check that we are in 2D
1266  IF(number_of_dimensions/=2) THEN
1267  local_error="The number of geometric dimensions of "// &
1268  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
1269  & " is invalid. The analytic function type of "// &
1270  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1271  & " requires that there be 2 geometric dimensions."
1272  CALL flagerror(local_error,err,error,*999)
1273  ENDIF
1274  !Create analytic field if required
1275  !Set analtyic function type
1276  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_hj_equation_two_dim_1
1277  CASE(equations_set_hj_equation_two_dim_2)
1278  !Check that we are in 2D
1279  IF(number_of_dimensions/=2) THEN
1280  local_error="The number of geometric dimensions of "// &
1281  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
1282  & " is invalid. The analytic function type of "// &
1283  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1284  & " requires that there be 2 geometric dimensions."
1285  CALL flagerror(local_error,err,error,*999)
1286  ENDIF
1287  !Create analytic field if required
1288  !Set analtyic function type
1289  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_hj_equation_two_dim_2
1290  CASE(equations_set_hj_equation_three_dim_1)
1291  !Check that we are in 3D
1292  IF(number_of_dimensions/=3) THEN
1293  local_error="The number of geometric dimensions of "// &
1294  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
1295  & " is invalid. The analytic function type of "// &
1296  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1297  & " requires that there be 3 geometric dimensions."
1298  CALL flagerror(local_error,err,error,*999)
1299  ENDIF
1300  !Create analytic field if required
1301  !Set analtyic function type
1302  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_hj_equation_three_dim_1
1303  CASE(equations_set_hj_equation_three_dim_2)
1304  !Check that we are in 3D
1305  IF(number_of_dimensions/=3) THEN
1306  local_error="The number of geometric dimensions of "// &
1307  & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
1308  & " is invalid. The analytic function type of "// &
1309  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1310  & " requires that there be 3 geometric dimensions."
1311  CALL flagerror(local_error,err,error,*999)
1312  ENDIF
1313  !Create analytic field if required
1314  !Set analtyic function type
1315  equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_hj_equation_three_dim_2
1316  CASE DEFAULT
1317  local_error="The specified analytic function type of "// &
1318  & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
1319  & " is invalid for a standard Hamilton-Jacobi equation."
1320  CALL flagerror(local_error,err,error,*999)
1321  END SELECT
1322  ELSE
1323  CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
1324  ENDIF
1325  ELSE
1326  CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
1327  ENDIF
1328  ELSE
1329  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
1330  ENDIF
1331  CASE(equations_set_setup_finish_action)
1332  IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
1333  analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
1334  IF(ASSOCIATED(analytic_field)) THEN
1335  IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
1336  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
1337  ENDIF
1338  ENDIF
1339  ELSE
1340  CALL flagerror("Equations set analytic is not associated.",err,error,*999)
1341  ENDIF
1342  CASE DEFAULT
1343  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1344  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1345  & " is invalid for a standard Hamilton-Jacobi equation."
1346  CALL flagerror(local_error,err,error,*999)
1347  END SELECT
1348  CASE(equations_set_setup_equations_type)
1349  SELECT CASE(equations_set_setup%ACTION_TYPE)
1350  CASE(equations_set_setup_start_action)
1351  IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
1352  CALL equations_create_start(equations_set,equations,err,error,*999)
1353  CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
1354  CALL equations_time_dependence_type_set(equations,equations_static,err,error,*999)
1355  ELSE
1356  CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
1357  ENDIF
1358  CASE(equations_set_setup_finish_action)
1359  SELECT CASE(equations_set%SOLUTION_METHOD)
1360  CASE(equations_set_fem_solution_method)
1361  !Finish the equations creation
1362  CALL equations_set_equations_get(equations_set,equations,err,error,*999)
1363  CALL equations_create_finish(equations,err,error,*999)
1364  !Create the equations mapping.
1365  CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
1366  CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
1367  CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,(/field_u_variable_type/), &
1368  & err,error,*999)
1369  CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
1370  CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
1371  !Create the equations matrices
1372  CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
1373  SELECT CASE(equations%SPARSITY_TYPE)
1374  CASE(equations_matrices_full_matrices)
1375  CALL equations_matrices_linear_storage_type_set(equations_matrices,(/matrix_block_storage_type/), &
1376  & err,error,*999)
1377  CASE(equations_matrices_sparse_matrices)
1378  CALL equations_matrices_linear_storage_type_set(equations_matrices,(/matrix_compressed_row_storage_type/), &
1379  & err,error,*999)
1380  CALL equationsmatrices_linearstructuretypeset(equations_matrices,(/equations_matrix_fem_structure/), &
1381  & err,error,*999)
1382  CASE DEFAULT
1383  local_error="The equations matrices sparsity type of "// &
1384  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
1385  CALL flagerror(local_error,err,error,*999)
1386  END SELECT
1387  CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
1388 ! CASE(EQUATIONS_SET_FMM_SOLUTION_METHOD)
1389 ! CALL FlagError("Not implemented.",ERR,ERROR,*999)
1390  CASE(equations_set_bem_solution_method)
1391  CALL flagerror("Not implemented.",err,error,*999)
1392  CASE(equations_set_fd_solution_method)
1393  CALL flagerror("Not implemented.",err,error,*999)
1394  CASE(equations_set_fv_solution_method)
1395  CALL flagerror("Not implemented.",err,error,*999)
1396  CASE(equations_set_gfem_solution_method)
1397  CALL flagerror("Not implemented.",err,error,*999)
1398  CASE(equations_set_gfv_solution_method)
1399  CALL flagerror("Not implemented.",err,error,*999)
1400  CASE DEFAULT
1401  local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
1402  & " is invalid."
1403  CALL flagerror(local_error,err,error,*999)
1404  END SELECT
1405  CASE DEFAULT
1406  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
1407  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1408  & " is invalid for a standard Hamilton-Jacobi equation."
1409  CALL flagerror(local_error,err,error,*999)
1410  END SELECT
1411  CASE DEFAULT
1412  local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
1413  & " is invalid for a standard Hamilton-Jacobi equation."
1414  CALL flagerror(local_error,err,error,*999)
1415  END SELECT
1416  ELSE
1417  local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
1418  & " does not equal a standard Hamilton-Jacobi equation subtype."
1419  CALL flagerror(local_error,err,error,*999)
1420  ENDIF
1421  ELSE
1422  CALL flagerror("Equations set is not associated.",err,error,*999)
1423  ENDIF
1424 
1425  exits("HJ_EQUATION_EQUATIONS_SET_STANDARD_SETUP")
1426  RETURN
1427 999 errorsexits("HJ_EQUATION_EQUATIONS_SET_STANDARD_SETUP",err,error)
1428  RETURN 1
1429  END SUBROUTINE hj_equation_equations_set_standard_setup
1430 
1431  !
1432  !================================================================================================================================
1433  !
1434 
1436  SUBROUTINE hj_equation_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
1437 
1438  !Argument variables
1439  TYPE(problem_type), POINTER :: PROBLEM
1440  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
1441  INTEGER(INTG), INTENT(OUT) :: ERR
1442  TYPE(varying_string), INTENT(OUT) :: ERROR
1443  !Local Variables
1444  TYPE(varying_string) :: LOCAL_ERROR
1445 
1446  enters("HJ_EQUATION_PROBLEM_SETUP",err,error,*999)
1447 
1448  IF(ASSOCIATED(problem)) THEN
1449  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
1450  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1451  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
1452  CALL flagerror("Problem specification must have three entries for a Hamilton-Jacobi problem.",err,error,*999)
1453  END IF
1454  SELECT CASE(problem%SPECIFICATION(3))
1455  CASE(problem_standard_hj_subtype)
1456  CALL hj_equation_problem_standard_setup(problem,problem_setup,err,error,*999)
1457  CASE(problem_generalised_hj_subtype)
1458  CALL flagerror("Not implemented.",err,error,*999)
1459  CASE DEFAULT
1460  local_error="Problem subtype "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
1461  & " is not valid for a Hamilton-Jacobi equation type of a classical field problem class."
1462  CALL flagerror(local_error,err,error,*999)
1463  END SELECT
1464  ELSE
1465  CALL flagerror("Problem is not associated.",err,error,*999)
1466  ENDIF
1467 
1468  exits("HJ_EQUATION_PROBLEM_SETUP")
1469  RETURN
1470 999 errorsexits("HJ_EQUATION_PROBLEM_SETUP",err,error)
1471  RETURN 1
1472  END SUBROUTINE hj_equation_problem_setup
1473 
1474  !
1475  !================================================================================================================================
1476  !
1477 
1479  SUBROUTINE hjequation_problemspecificationset(problem,problemSpecification,err,error,*)
1480 
1481  !Argument variables
1482  TYPE(problem_type), POINTER :: problem
1483  INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
1484  INTEGER(INTG), INTENT(OUT) :: err
1485  TYPE(varying_string), INTENT(OUT) :: error
1486  !Local Variables
1487  TYPE(varying_string) :: localError
1488  INTEGER(INTG) :: problemSubtype
1489 
1490  enters("HJEquation_ProblemSpecificationSet",err,error,*999)
1491 
1492  IF(ASSOCIATED(problem)) THEN
1493  IF(SIZE(problemspecification,1)==3) THEN
1494  problemsubtype=problemspecification(3)
1495  SELECT CASE(problemsubtype)
1496  CASE(problem_standard_hj_subtype)
1497  !ok
1498  CASE(problem_generalised_hj_subtype)
1499  CALL flagerror("Not implemented.",err,error,*999)
1500  CASE DEFAULT
1501  localerror="The third problem specification of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
1502  & " is not valid for a Hamilton-Jacobi type of a classical field problem."
1503  CALL flagerror(localerror,err,error,*999)
1504  END SELECT
1505  IF(ALLOCATED(problem%specification)) THEN
1506  CALL flagerror("Problem specification is already allocated.",err,error,*999)
1507  ELSE
1508  ALLOCATE(problem%specification(3),stat=err)
1509  IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
1510  END IF
1511  problem%specification(1:3)=[problem_classical_field_class,problem_hj_equation_type,problemsubtype]
1512  ELSE
1513  CALL flagerror("Hamilton-Jacobi problem specification must have three entries.",err,error,*999)
1514  END IF
1515  ELSE
1516  CALL flagerror("Problem is not associated.",err,error,*999)
1517  END IF
1518 
1519  exits("HJEquation_ProblemSpecificationSet")
1520  RETURN
1521 999 errors("HJEquation_ProblemSpecificationSet",err,error)
1522  exits("HJEquation_ProblemSpecificationSet")
1523  RETURN 1
1524 
1525  END SUBROUTINE hjequation_problemspecificationset
1526 
1527  !
1528  !================================================================================================================================
1529  !
1530 
1532  SUBROUTINE hj_equation_problem_standard_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
1533 
1534  !Argument variables
1535  TYPE(problem_type), POINTER :: PROBLEM
1536  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
1537  INTEGER(INTG), INTENT(OUT) :: ERR
1538  TYPE(varying_string), INTENT(OUT) :: ERROR
1539  !Local Variables
1540  TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
1541  TYPE(solver_type), POINTER :: SOLVER
1542  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
1543  TYPE(solvers_type), POINTER :: SOLVERS
1544  TYPE(varying_string) :: LOCAL_ERROR
1545 
1546  enters("HJ_EQUATION_PROBLEM_STANDARD_SETUP",err,error,*999)
1547 
1548  NULLIFY(control_loop)
1549  NULLIFY(solver)
1550  NULLIFY(solver_equations)
1551  NULLIFY(solvers)
1552  IF(ASSOCIATED(problem)) THEN
1553  IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
1554  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1555  ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
1556  CALL flagerror("Problem specification must have three entries for a Hamilton-Jacobi problem.",err,error,*999)
1557  END IF
1558  IF(problem%SPECIFICATION(3)==problem_standard_hj_subtype) THEN
1559  SELECT CASE(problem_setup%SETUP_TYPE)
1560  CASE(problem_setup_initial_type)
1561  SELECT CASE(problem_setup%ACTION_TYPE)
1562  CASE(problem_setup_start_action)
1563  !Do nothing????
1564  CASE(problem_setup_finish_action)
1565  !Do nothing???
1566  CASE DEFAULT
1567  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1568  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1569  & " is invalid for a standard Hamilton-Jacobi equation."
1570  CALL flagerror(local_error,err,error,*999)
1571  END SELECT
1572  CASE(problem_setup_control_type)
1573  SELECT CASE(problem_setup%ACTION_TYPE)
1574  CASE(problem_setup_start_action)
1575  !Set up a simple control loop
1576  CALL control_loop_create_start(problem,control_loop,err,error,*999)
1577  CASE(problem_setup_finish_action)
1578  !Finish the control loops
1579  control_loop_root=>problem%CONTROL_LOOP
1580  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1581  CALL control_loop_create_finish(control_loop,err,error,*999)
1582  CASE DEFAULT
1583  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1584  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1585  & " is invalid for a standard Hamilton-Jacobi equation."
1586  CALL flagerror(local_error,err,error,*999)
1587  END SELECT
1588  CASE(problem_setup_solvers_type)
1589  !Get the control loop
1590  control_loop_root=>problem%CONTROL_LOOP
1591  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1592  SELECT CASE(problem_setup%ACTION_TYPE)
1593  CASE(problem_setup_start_action)
1594  !Start the solvers creation
1595  CALL solvers_create_start(control_loop,solvers,err,error,*999)
1596  CALL solvers_number_set(solvers,1,err,error,*999)
1597  !Set the solver to be a linear solver
1598  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1599  CALL solver_type_set(solver,solver_linear_type,err,error,*999)
1600  !Set solver defaults
1601  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
1602  CASE(problem_setup_finish_action)
1603  !Get the solvers
1604  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1605  !Finish the solvers creation
1606  CALL solvers_create_finish(solvers,err,error,*999)
1607  CASE DEFAULT
1608  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1609  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1610  & " is invalid for a standard Hamilton-Jacobi equation."
1611  CALL flagerror(local_error,err,error,*999)
1612  END SELECT
1613  CASE(problem_setup_solver_equations_type)
1614  SELECT CASE(problem_setup%ACTION_TYPE)
1615  CASE(problem_setup_start_action)
1616  !Get the control loop
1617  control_loop_root=>problem%CONTROL_LOOP
1618  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1619  !Get the solver
1620  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1621  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1622  !Create the solver equations
1623  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
1624  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
1625  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_static,err,error,*999)
1626  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
1627  CASE(problem_setup_finish_action)
1628  !Get the control loop
1629  control_loop_root=>problem%CONTROL_LOOP
1630  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
1631  !Get the solver equations
1632  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1633  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1634  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1635  !Finish the solver equations creation
1636  CALL solver_equations_create_finish(solver_equations,err,error,*999)
1637  CASE DEFAULT
1638  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
1639  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1640  & " is invalid for a standard Hamilton-Jacobi equation."
1641  CALL flagerror(local_error,err,error,*999)
1642  END SELECT
1643  CASE DEFAULT
1644  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
1645  & " is invalid for a standard Hamilton-Jacobi equation."
1646  CALL flagerror(local_error,err,error,*999)
1647  END SELECT
1648  ELSE
1649  local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
1650  & " does not equal a standard Hamilton-Jacobi equation subtype."
1651  CALL flagerror(local_error,err,error,*999)
1652  ENDIF
1653  ELSE
1654  CALL flagerror("Problem is not associated.",err,error,*999)
1655  ENDIF
1656 
1657  exits("HJ_EQUATION_PROBLEM_STANDARD_SETUP")
1658  RETURN
1659 999 errorsexits("HJ_EQUATION_PROBLEM_STANDARD_SETUP",err,error)
1660  RETURN 1
1661  END SUBROUTINE hj_equation_problem_standard_setup
1662 
1663  !
1664  !================================================================================================================================
1665  !
1666 !
1667 !
1668 !
1669 
1670 
1672  SUBROUTINE number_of_input_nodes(INPUT_FILE_NAME,INPUT_FILE_FORMAT,TOTAL_NUMBER_OF_NODES,TOTAL_NUMBER_OF_ELEMENTS,&
1673  &total_number_of_connectivity,err)
1674 
1675  !subroutine variables
1676  INTEGER(INTG), INTENT(OUT) :: TOTAL_NUMBER_OF_NODES
1677  INTEGER(INTG), INTENT(OUT) :: TOTAL_NUMBER_OF_ELEMENTS
1678  INTEGER(INTG), INTENT(OUT) :: TOTAL_NUMBER_OF_CONNECTIVITY
1679  CHARACTER (LEN=300) :: INPUT_FILE_NAME
1680  CHARACTER (LEN=10) :: INPUT_FILE_FORMAT
1681  INTEGER(INTG) :: Err
1682 
1683  !Argument variables
1684 ! TYPE(VARYING_STRING) :: LOCAL_ERROR !<The error string
1685 
1686  !Local variables
1687  INTEGER(INTG), ALLOCATABLE, DIMENSION(:,:):: CONNECTIVITY_LIST
1688  INTEGER(INTG), ALLOCATABLE, DIMENSION(:,:):: ELEMENT_LIST
1689  INTEGER(INTG), ALLOCATABLE, DIMENSION(:) :: CONNECTIVITY_NUMBER
1690  INTEGER(INTG) :: I,J,K,N,NUMBER_OF_NODES_PER_ELEMENT,THERE_IS_IN_CONNECTIVITY_LIST
1691  CHARACTER (LEN=300) :: STRING
1692 
1693 ! ENTERS("GENERATE_STATUS_MASK",Err,Error,*999)
1694 
1695 
1696 ! """""""""""""""""""""""""""""""""""INPUT OF TABC FORMAT""""""""""""""""""""""""""""""
1697  IF (input_file_format .EQ. "TABC") THEN
1698 
1699  i = index(input_file_name,' ') - 1
1700  string = input_file_name(1:i)//".tabc"
1701  OPEN (11,file=string)
1702  READ(11,*) string
1703 ! PRINT *, STRING
1704  total_number_of_nodes=-1
1705  DO WHILE (string .ne. "Connectivity")
1706  READ(11,*) string
1707  total_number_of_nodes=total_number_of_nodes+1
1708  ENDDO
1709 
1710  total_number_of_connectivity=0
1711  DO i=1,total_number_of_nodes
1712  READ(11,*) string,j
1713  total_number_of_connectivity=total_number_of_connectivity+j
1714  ENDDO
1715 
1716  CLOSE (11)
1717  total_number_of_elements = total_number_of_connectivity ! SHOULD BE DEFINED
1718  ENDIF
1719 
1720 ! """""""""""""""""""""""""""""""""""INPUT OF VTK FORMAT""""""""""""""""""""""""""""""
1721  IF (input_file_format .EQ. "VTKTET") THEN
1722 
1723  number_of_nodes_per_element = 4
1724 
1725  i = index(input_file_name,' ') - 1
1726  string = input_file_name(1:i)//".vtk"
1727 ! PRINT *, STRING
1728  OPEN (11,file=string)
1729  READ(11,*) string
1730  READ(11,*) string
1731  READ(11,*) string
1732  READ(11,*) string
1733  READ(11,*) string,total_number_of_nodes
1734 ! PRINT *, STRING,TOTAL_NUMBER_OF_NODES
1735 
1736  DO i=1,int((total_number_of_nodes/3.0_dp)+0.5_dp)
1737 
1738  READ(11,*) string
1739 
1740  ENDDO
1741 ! READ(11,*) STRING
1742 ! print*,I,INT(TOTAL_NUMBER_OF_NODES/3.0+0.5),STRING
1743  READ(11,*) string,total_number_of_elements
1744 
1745  ALLOCATE(element_list(total_number_of_elements,20),stat=err)
1746  ALLOCATE(connectivity_list(total_number_of_nodes,50),stat=err)
1747  ALLOCATE(connectivity_number(total_number_of_nodes),stat=err)
1748 
1749  DO i=1,total_number_of_elements
1750  DO j=1,20
1751  element_list(i,j)=0
1752  ENDDO
1753  ENDDO
1754  DO i=1,total_number_of_nodes
1755  connectivity_number(i)=0
1756  DO j=1,50
1757  connectivity_list(i,j)=0
1758  ENDDO
1759  ENDDO
1760 
1761  total_number_of_connectivity=0
1762 
1763  DO i=1,total_number_of_elements
1764 
1765  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
1766 
1767 ! we have nodes starting number of 0
1768  DO j=1,number_of_nodes_per_element
1769  element_list(i,j)=element_list(i,j)+1
1770  ENDDO
1771 
1772  DO j=1,number_of_nodes_per_element-1
1773 
1774  DO k=1+j,number_of_nodes_per_element
1775  there_is_in_connectivity_list = 0
1776  DO n=1,connectivity_number(element_list(i,j))
1777 
1778  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
1779  there_is_in_connectivity_list = 1
1780  ENDIF
1781 
1782  ENDDO
1783 
1784  IF (there_is_in_connectivity_list .EQ. 0) THEN
1785 
1786  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
1787  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
1788  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
1789  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
1790 
1791  total_number_of_connectivity=total_number_of_connectivity+1
1792 
1793  ENDIF
1794 
1795  ENDDO
1796 
1797  ENDDO
1798 
1799  ENDDO
1800  total_number_of_connectivity=total_number_of_connectivity*2
1801 
1802  CLOSE (11)
1803 
1804  ENDIF
1805 ! """""""""""""""""""""""""""""""""""INPUT OF VTK FORMAT""""""""""""""""""""""""""""""
1806  IF (input_file_format .EQ. "VTKTET1NPL") THEN
1807 
1808  number_of_nodes_per_element = 4
1809 
1810  i = index(input_file_name,' ') - 1
1811  string = input_file_name(1:i)//".vtk"
1812 ! PRINT *, STRING
1813  OPEN (11,file=string)
1814  READ(11,*) string
1815  READ(11,*) string
1816  READ(11,*) string
1817  READ(11,*) string
1818  READ(11,*) string,total_number_of_nodes,string
1819 ! PRINT *, STRING,TOTAL_NUMBER_OF_NODES,STRING
1820 
1821  DO i=1,total_number_of_nodes
1822 
1823  READ(11,*) string
1824 
1825  ENDDO
1826 ! READ(11,*) STRING
1827  READ(11,*) string,total_number_of_elements
1828 
1829  ALLOCATE(element_list(total_number_of_elements,20),stat=err)
1830  ALLOCATE(connectivity_list(total_number_of_nodes,50),stat=err)
1831  ALLOCATE(connectivity_number(total_number_of_nodes),stat=err)
1832 
1833  DO i=1,total_number_of_elements
1834  DO j=1,20
1835  element_list(i,j)=0
1836  ENDDO
1837  ENDDO
1838  DO i=1,total_number_of_nodes
1839  connectivity_number(i)=0
1840  DO j=1,50
1841  connectivity_list(i,j)=0
1842  ENDDO
1843  ENDDO
1844 
1845  total_number_of_connectivity=0
1846 
1847  DO i=1,total_number_of_elements
1848 
1849  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
1850 
1851 ! we have nodes starting number of 0
1852  DO j=1,number_of_nodes_per_element
1853  element_list(i,j)=element_list(i,j)+1
1854  ENDDO
1855 
1856  DO j=1,number_of_nodes_per_element-1
1857 
1858  DO k=1+j,number_of_nodes_per_element
1859  there_is_in_connectivity_list = 0
1860  DO n=1,connectivity_number(element_list(i,j))
1861 
1862  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
1863  there_is_in_connectivity_list = 1
1864  ENDIF
1865 
1866  ENDDO
1867 
1868  IF (there_is_in_connectivity_list .EQ. 0) THEN
1869 
1870  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
1871  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
1872  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
1873  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
1874 
1875  total_number_of_connectivity=total_number_of_connectivity+1
1876 
1877  ENDIF
1878 
1879  ENDDO
1880 
1881  ENDDO
1882 
1883  ENDDO
1884  total_number_of_connectivity=total_number_of_connectivity*2
1885 
1886  CLOSE (11)
1887 
1888  ENDIF
1889 ! """""""""""""""""""""""""""""""""""INPUT OF CARP FORMAT""""""""""""""""""""""""""""""
1890  IF (input_file_format .EQ. "CARP") THEN
1891 
1892  number_of_nodes_per_element = 4
1893 
1894  i = index(input_file_name,' ') - 1
1895  string = input_file_name(1:i)//".pts"
1896 ! PRINT *, STRING
1897  OPEN (11,file=string)
1898  READ(11,*) total_number_of_nodes
1899 ! PRINT *, TOTAL_NUMBER_OF_NODES
1900  CLOSE (11)
1901 
1902  string = input_file_name(1:i)//".elem"
1903 ! PRINT *, STRING
1904  OPEN (11,file=string)
1905  READ(11,*) total_number_of_elements
1906 
1907  ALLOCATE(element_list(total_number_of_elements,20),stat=err)
1908  ALLOCATE(connectivity_list(total_number_of_nodes,50),stat=err)
1909  ALLOCATE(connectivity_number(total_number_of_nodes),stat=err)
1910 
1911  DO i=1,total_number_of_elements
1912  DO j=1,20
1913  element_list(i,j)=0
1914  ENDDO
1915  ENDDO
1916  DO i=1,total_number_of_nodes
1917  connectivity_number(i)=0
1918  DO j=1,50
1919  connectivity_list(i,j)=0
1920  ENDDO
1921  ENDDO
1922 
1923  total_number_of_connectivity=0
1924 
1925  DO i=1,total_number_of_elements
1926 
1927  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
1928 
1929 ! we have nodes starting number of 0
1930  DO j=1,number_of_nodes_per_element
1931  element_list(i,j)=element_list(i,j)+1
1932  ENDDO
1933 
1934  DO j=1,number_of_nodes_per_element-1
1935 
1936  DO k=1+j,number_of_nodes_per_element
1937  there_is_in_connectivity_list = 0
1938  DO n=1,connectivity_number(element_list(i,j))
1939 
1940  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
1941  there_is_in_connectivity_list = 1
1942  ENDIF
1943 
1944  ENDDO
1945 
1946  IF (there_is_in_connectivity_list .EQ. 0) THEN
1947 
1948  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
1949  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
1950  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
1951  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
1952 
1953  total_number_of_connectivity=total_number_of_connectivity+1
1954 
1955  ENDIF
1956 
1957  ENDDO
1958 
1959  ENDDO
1960 
1961  ENDDO
1962  total_number_of_connectivity=total_number_of_connectivity*2
1963 
1964  CLOSE (11)
1965 
1966  ENDIF
1967 
1968 ! """""""""""""""""""""""""""""""""""INPUT OF TETGEN FORMAT""""""""""""""""""""""""""""""
1969  IF (input_file_format .EQ. "TETGEN") THEN
1970 
1971  number_of_nodes_per_element = 4
1972 
1973  i = index(input_file_name,' ') - 1
1974  string = input_file_name(1:i)//".node"
1975  OPEN (11,file=string)
1976  READ(11,*) total_number_of_nodes
1977  CLOSE (11)
1978 
1979  string = input_file_name(1:i)//".ele"
1980  OPEN (11,file=string)
1981  READ(11,*) total_number_of_elements
1982 
1983  ALLOCATE(element_list(total_number_of_elements,20),stat=err)
1984  ALLOCATE(connectivity_list(total_number_of_nodes,50),stat=err)
1985  ALLOCATE(connectivity_number(total_number_of_nodes),stat=err)
1986 
1987  DO i=1,total_number_of_elements
1988  DO j=1,20
1989  element_list(i,j)=0
1990  ENDDO
1991  ENDDO
1992  DO i=1,total_number_of_nodes
1993  connectivity_number(i)=0
1994  DO j=1,50
1995  connectivity_list(i,j)=0
1996  ENDDO
1997  ENDDO
1998 
1999  total_number_of_connectivity=0
2000 
2001  DO i=1,total_number_of_elements
2002 
2003  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
2004 
2005  DO j=1,number_of_nodes_per_element-1
2006 
2007  DO k=1+j,number_of_nodes_per_element
2008  there_is_in_connectivity_list = 0
2009  DO n=1,connectivity_number(element_list(i,j))
2010 
2011  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
2012  there_is_in_connectivity_list = 1
2013  ENDIF
2014 
2015  ENDDO
2016 
2017  IF (there_is_in_connectivity_list .EQ. 0) THEN
2018  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
2019  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
2020  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
2021  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
2022 
2023  total_number_of_connectivity=total_number_of_connectivity+1
2024 
2025  ENDIF
2026 
2027  ENDDO
2028 
2029  ENDDO
2030 
2031  ENDDO
2032 
2033  CLOSE (11)
2034 
2035  total_number_of_connectivity=total_number_of_connectivity*2
2036 
2037  ENDIF
2038 
2039 
2040  END SUBROUTINE number_of_input_nodes
2041 
2042  !
2043  !================================================================================================================================
2044  !
2045 
2046 
2048  SUBROUTINE pre_process_information(MATERIAL_BEHAVIOUR,INPUT_FILE_NAME,INPUT_FILE_FORMAT,TOTAL_NUMBER_OF_NODES,&
2049 &input_type_for_seed_value,input_type_for_speed_function,speed_function_along_eigen_vector,input_type_for_conductivity,&
2050 &status_mask,node_list,conductivity_tensor,speed_function_table,seed_value,connectivity_number,&
2051 &speed_function_table_on_connectivity,conductivity_tensor_on_connectivity,raw_index,column_index,total_number_of_connectivity,&
2052 &connectivity_list,element_list,total_number_of_elements,number_of_nodes_per_element,err)
2053 
2054  !subroutine variables
2055  REAL(DP), ALLOCATABLE :: NODE_LIST(:,:)
2056  REAL(DP), ALLOCATABLE :: SPEED_FUNCTION_TABLE(:,:)
2057  REAL(DP), ALLOCATABLE :: CONDUCTIVITY_TENSOR(:,:)
2058  REAL(DP), ALLOCATABLE :: SPEED_FUNCTION_TABLE_ON_CONNECTIVITY(:,:)
2059  REAL(DP), ALLOCATABLE :: CONDUCTIVITY_TENSOR_ON_CONNECTIVITY(:,:)
2060  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_LIST(:,:)
2061  INTEGER(INTG), ALLOCATABLE :: ELEMENT_LIST(:,:)
2062 
2063  INTEGER(INTG), ALLOCATABLE, DIMENSION(:) :: COLUMN_INDEX
2064  INTEGER(INTG), ALLOCATABLE, DIMENSION(:) :: RAW_INDEX
2065 
2066  REAL(DP), ALLOCATABLE :: SEED_VALUE(:)
2067  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_NUMBER(:)
2068  CHARACTER (LEN=10), ALLOCATABLE :: STATUS_MASK(:)
2069  REAL(DP) :: SPEED_FUNCTION_ALONG_EIGEN_VECTOR(3)
2070  INTEGER(INTG), INTENT(OUT) :: TOTAL_NUMBER_OF_ELEMENTS,TOTAL_NUMBER_OF_CONNECTIVITY
2071  INTEGER(INTG), INTENT(OUT) :: NUMBER_OF_NODES_PER_ELEMENT
2072  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_NODES
2073  CHARACTER (LEN=10) :: INPUT_TYPE_FOR_SEED_VALUE
2074  CHARACTER (LEN=10) :: INPUT_TYPE_FOR_SPEED_FUNCTION
2075  CHARACTER (LEN=10) :: INPUT_TYPE_FOR_CONDUCTIVITY
2076  CHARACTER (LEN=300) :: INPUT_FILE_NAME
2077  CHARACTER (LEN=10) :: INPUT_FILE_FORMAT
2078  CHARACTER (LEN=12) :: MATERIAL_BEHAVIOUR
2079  INTEGER(INTG) :: Err
2080  TYPE(varying_string) :: Error
2081 
2082  !Local variables
2083  CHARACTER (LEN=300) :: STRING
2084  INTEGER(INTG) :: I,J,K,N,FIRST_NODE_NUMBER
2085  INTEGER(INTG) :: TEXT_LENGTH, THERE_IS_IN_CONNECTIVITY_LIST
2086  REAL(DP) :: A(3),B(3),C(3)
2087  REAL(DP) :: DOT_PRODUCT_VALUE
2088 
2089 !INITIALIZE PARAMETERS:
2090  DO i=1,total_number_of_nodes
2091  connectivity_number(i)=0
2092  DO j=1,3
2093  node_list(i,j) = 0.0
2094  speed_function_table(i,j) = 0.0
2095  ENDDO
2096  DO j=1,9
2097  conductivity_tensor(i,j) = 0.0
2098  ENDDO
2099  raw_index(i)=0
2100  ENDDO
2101  raw_index(total_number_of_nodes+1)=0
2102 
2103  DO i=1,total_number_of_connectivity
2104  column_index(i)=0
2105  DO j=1,3
2106  speed_function_table_on_connectivity(i,j) = 0.0
2107  ENDDO
2108  DO j=1,9
2109  conductivity_tensor_on_connectivity(i,j) = 0.0
2110  ENDDO
2111  ENDDO
2112 
2113 ! """""""""""""""""""""""""""""""""""INPUT OF TABC FORMAT""""""""""""""""""""""""""""""
2114  IF (input_file_format .EQ. "TABC") THEN
2115 
2116  i = index(input_file_name,' ') - 1
2117  string = input_file_name(1:i)//".tabc"
2118  OPEN (11,file=string)
2119  READ(11,*) string
2120 
2121 ! SOSIOISOISOISOISOIS load data for the case material behaves = ISOTROPIC
2122  IF (material_behaviour .EQ. "ISOTROPIC") THEN
2123 
2124 ! input type for *velocity function* = FILE and *seed points* = FILE
2125  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "FILE") THEN
2126  DO i=1,total_number_of_nodes
2127 
2128  READ(11,*) string,(node_list(i,j),j=1,3),speed_function_along_eigen_vector(1),seed_value(i)
2129 
2130  speed_function_along_eigen_vector(2)=speed_function_along_eigen_vector(1)
2131  speed_function_along_eigen_vector(3)=speed_function_along_eigen_vector(1)
2132 
2133  conductivity_tensor(i,1)=1.0_dp
2134  conductivity_tensor(i,2)=0.0_dp
2135  conductivity_tensor(i,3)=0.0_dp
2136  conductivity_tensor(i,4)=0.0_dp
2137  conductivity_tensor(i,5)=1.0_dp
2138  conductivity_tensor(i,6)=0.0_dp
2139  conductivity_tensor(i,7)=0.0_dp
2140  conductivity_tensor(i,8)=0.0_dp
2141  conductivity_tensor(i,9)=1.0_dp
2142 
2143  ENDDO
2144  ENDIF
2145 
2146 ! input type for *velocity function* = FIXED and *seed points* = FILE
2147  IF (input_type_for_speed_function .EQ. "FIXED" .AND. input_type_for_seed_value .EQ. "FILE") THEN
2148  DO i=1,total_number_of_nodes
2149 
2150  READ(11,*) string,(node_list(i,j),j=1,3),seed_value(i)
2151  conductivity_tensor(i,1)=1.0_dp
2152  conductivity_tensor(i,2)=0.0_dp
2153  conductivity_tensor(i,3)=0.0_dp
2154  conductivity_tensor(i,4)=0.0_dp
2155  conductivity_tensor(i,5)=1.0_dp
2156  conductivity_tensor(i,6)=0.0_dp
2157  conductivity_tensor(i,7)=0.0_dp
2158  conductivity_tensor(i,8)=0.0_dp
2159  conductivity_tensor(i,9)=1.0_dp
2160 
2161  speed_function_along_eigen_vector(2)=speed_function_along_eigen_vector(1)
2162  speed_function_along_eigen_vector(3)=speed_function_along_eigen_vector(1)
2163 
2164  ENDDO
2165  ENDIF
2166 
2167 ! input type for *velocity function* = FILE and *seed points* = LIST
2168  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "LIST") THEN
2169  DO i=1,total_number_of_nodes
2170 
2171  READ(11,*) string,(node_list(i,j),j=1,3),speed_function_along_eigen_vector(1)
2172  conductivity_tensor(i,1)=1.0_dp
2173  conductivity_tensor(i,2)=0.0_dp
2174  conductivity_tensor(i,3)=0.0_dp
2175  conductivity_tensor(i,4)=0.0_dp
2176  conductivity_tensor(i,5)=1.0_dp
2177  conductivity_tensor(i,6)=0.0_dp
2178  conductivity_tensor(i,7)=0.0_dp
2179  conductivity_tensor(i,8)=0.0_dp
2180  conductivity_tensor(i,9)=1.0_dp
2181 
2182  speed_function_along_eigen_vector(2)=speed_function_along_eigen_vector(1)
2183  speed_function_along_eigen_vector(3)=speed_function_along_eigen_vector(1)
2184 
2185  IF (status_mask(i) .NE. "SEED POINT") THEN
2186  seed_value(i) = 1000.0_dp
2187  ENDIF
2188 
2189  ENDDO
2190  ENDIF
2191 
2192 ! input type for *velocity function* = FIXED and *seed points* = LIST
2193  IF (input_type_for_speed_function .EQ. "FIXED" .AND. input_type_for_seed_value .EQ. "LIST") THEN
2194  DO i=1,total_number_of_nodes
2195 
2196  READ(11,*) string,(node_list(i,j),j=1,3)
2197  conductivity_tensor(i,1)=1.0_dp
2198  conductivity_tensor(i,2)=0.0_dp
2199  conductivity_tensor(i,3)=0.0_dp
2200  conductivity_tensor(i,4)=0.0_dp
2201  conductivity_tensor(i,5)=1.0_dp
2202  conductivity_tensor(i,6)=0.0_dp
2203  conductivity_tensor(i,7)=0.0_dp
2204  conductivity_tensor(i,8)=0.0_dp
2205  conductivity_tensor(i,9)=1.0_dp
2206 
2207  speed_function_along_eigen_vector(2)=speed_function_along_eigen_vector(1)
2208  speed_function_along_eigen_vector(3)=speed_function_along_eigen_vector(1)
2209 
2210  IF (status_mask(i) .NE. "SEED POINT") THEN
2211  seed_value(i) = 1000.0_dp
2212  ENDIF
2213 
2214  ENDDO
2215  ENDIF
2216 
2217  ENDIF
2218 
2219 ! ANANANAIANSOANAIANI load data for the case material behaves = ANISOTROPIC
2220  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
2221 
2222 
2223 ! if conductivity format is TENSOR type i.e. three EigenVectors
2224  IF (input_type_for_conductivity .EQ. "TENSOR") THEN
2225 
2226 ! input type for *velocity function* = FILE and *seed points* = FILE
2227  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "FILE") THEN
2228  DO i=1,total_number_of_nodes
2229 
2230  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,9),speed_function_along_eigen_vector(1),&
2231  &speed_function_along_eigen_vector(2),speed_function_along_eigen_vector(3),seed_value(i)
2232 
2233  ENDDO
2234  ENDIF
2235 
2236 ! input type for *velocity function* = FIXED and *seed points* = FILE
2237  IF (input_type_for_speed_function .EQ. "FIXED" .AND. input_type_for_seed_value .EQ. "FILE") THEN
2238  DO i=1,total_number_of_nodes
2239 
2240  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,9),seed_value(i)
2241 
2242  ENDDO
2243  ENDIF
2244 
2245 ! input type for *velocity function* = FILE and *seed points* = LIST
2246  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "LIST") THEN
2247  DO i=1,total_number_of_nodes
2248 
2249  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,9),speed_function_along_eigen_vector(1),&
2250  &speed_function_along_eigen_vector(2),speed_function_along_eigen_vector(3)
2251 
2252  IF (status_mask(i) .NE. "SEED POINT") THEN
2253  seed_value(i) = 1000.0_dp
2254  ENDIF
2255 
2256  ENDDO
2257  ENDIF
2258 
2259 ! input type for *velocity function* = FIXED and *seed points* = LIST
2260  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "LIST") THEN
2261  DO i=1,total_number_of_nodes
2262 
2263  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,9)
2264 
2265  IF (status_mask(i) .NE. "SEED POINT") THEN
2266  seed_value(i) = 1000.0_dp
2267  ENDIF
2268 
2269  ENDDO
2270  ENDIF
2271 
2272  ENDIF
2273 
2274 
2275 ! if conductivity format is VECTOR type i.e. first EigenVectors
2276  IF (input_type_for_conductivity .EQ. "VECTOR") THEN
2277 
2278 ! input type for *velocity function* = FILE and *seed points* = FILE
2279  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "FILE") THEN
2280  DO i=1,total_number_of_nodes
2281 
2282  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,3),speed_function_along_eigen_vector(1),&
2283  &speed_function_along_eigen_vector(2),speed_function_along_eigen_vector(3),seed_value(i)
2284 
2285  ENDDO
2286  ENDIF
2287 
2288 ! input type for *velocity function* = FIXED and *seed points* = FILE
2289  IF (input_type_for_speed_function .EQ. "FIXED" .AND. input_type_for_seed_value .EQ. "FILE") THEN
2290  DO i=1,total_number_of_nodes
2291 
2292  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,3),seed_value(i)
2293 
2294  ENDDO
2295  ENDIF
2296 
2297 ! input type for *velocity function* = FILE and *seed points* = LIST
2298  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "LIST") THEN
2299  DO i=1,total_number_of_nodes
2300 
2301  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,3),speed_function_along_eigen_vector(1),&
2302  &speed_function_along_eigen_vector(2),speed_function_along_eigen_vector(3)
2303 
2304  IF (status_mask(i) .NE. "SEED POINT") THEN
2305  seed_value(i) = 1000.0_dp
2306  ENDIF
2307 
2308  ENDDO
2309  ENDIF
2310 
2311 ! input type for *velocity function* = FIXED and *seed points* = LIST
2312  IF (input_type_for_speed_function .EQ. "FILE" .AND. input_type_for_seed_value .EQ. "LIST") THEN
2313  DO i=1,total_number_of_nodes
2314 
2315  READ(11,*) string,(node_list(i,j),j=1,3),(conductivity_tensor(i,j),j=1,3)
2316 
2317  IF (status_mask(i) .NE. "SEED POINT") THEN
2318  seed_value(i) = 1000.0_dp
2319  ENDIF
2320 
2321  ENDDO
2322  ENDIF
2323 
2324  DO i=1,total_number_of_nodes
2325 ! CALL CALCULATE_SECOND_EIGENVECTOR()
2326  a=(/conductivity_tensor(i,1),conductivity_tensor(i,2),conductivity_tensor(i,3)/)
2327  b=(/0.0_dp,0.0_dp,1.0_dp/)
2328  CALL vector_vector_product(a,b,dot_product_value,err)
2329 
2330  IF (dot_product_value .LT. 0.8_dp) THEN
2331  CALL cross_product(a,b,c,err,error,*999)
2332  ELSE
2333  b=(/0.0_dp,1.0_dp,0.0_dp/)
2334  CALL vector_vector_product(a,b,dot_product_value,err)
2335  IF (dot_product_value .LT. 0.8_dp) THEN
2336  CALL cross_product(a,b,c,err,error,*999)
2337  ELSE
2338  b=(/1.0_dp,0.0_dp,0.0_dp/)
2339  CALL cross_product(a,b,c,err,error,*999)
2340  ENDIF
2341  ENDIF
2342 
2343  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>=zero_tolerance) THEN
2344 
2345  conductivity_tensor(i,4) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2346  conductivity_tensor(i,5) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2347  conductivity_tensor(i,6) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2348 
2349  ENDIF
2350 
2351  b=(/conductivity_tensor(i,4),conductivity_tensor(i,5),conductivity_tensor(i,6)/)
2352  CALL cross_product(a,b,c,err,error,*999)
2353 
2354  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>zero_tolerance) THEN
2355 
2356  conductivity_tensor(i,7) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2357  conductivity_tensor(i,8) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2358  conductivity_tensor(i,9) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2359 
2360  ENDIF
2361 
2362 ! CALL CALCULATE_SECOND_EIGENVECTOR()
2363  ENDDO
2364 
2365  ENDIF
2366 
2367  ENDIF
2368 
2369 ! CONNSDONCOCNCNOSKCN load data for the CONNECTIVITY list
2370  READ(11,*) string
2371 
2372  DO i=1,total_number_of_nodes
2373 
2374  READ(11,*) string,connectivity_number(i),(connectivity_list(i,j),j=1,connectivity_number(i))
2375 
2376  DO j=1,3
2377  speed_function_table(i,j)=speed_function_along_eigen_vector(j)
2378  ENDDO
2379 
2380 ! PRINT *,STRING,CONNECTIVITY_NUMBER(I),(CONNECTIVITY_LIST(I,J),J=1,CONNECTIVITY_NUMBER(I))
2381 
2382  ENDDO
2383 
2384  CLOSE(11)
2385 
2386  ENDIF
2387 
2388 
2389 ! """""""""""""""""""""""""""""""""""INPUT OF VTK FORMAT""""""""""""""""""""""""""""""
2390  IF (input_file_format .EQ. "VTKTET") THEN
2391 
2392  number_of_nodes_per_element = 4
2393 
2394 ! NDNAPDNDONOEEENODED load data for NODAL POSITIONING info
2395  text_length = index(input_file_name,' ') - 1
2396  string = input_file_name(1:text_length)//".vtk"
2397 
2398  OPEN (11,file=string)
2399  READ(11,*) string
2400  READ(11,*) string
2401  READ(11,*) string
2402  READ(11,*) string
2403  READ(11,*) string
2404 
2405  DO i=1,int((total_number_of_nodes/3.0_dp)+0.5_dp)-1
2406 
2407  READ(11,*) (node_list(3*(i-1)+1,j),j=1,3),(node_list(3*(i-1)+2,j),j=1,3),(node_list(3*(i-1)+3,j),j=1,3)
2408 
2409  ENDDO
2410 
2411  i=int((total_number_of_nodes/3.0_dp)+0.5_dp)
2412 
2413  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 0) THEN
2414 
2415  READ(11,*) (node_list(3*(i-1)+1,j),j=1,3),(node_list(3*(i-1)+2,j),j=1,3),(node_list(3*(i-1)+3,j),j=1,3)
2416 
2417  ENDIF
2418  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 1) THEN
2419 
2420  READ(11,*) (node_list(3*(i-1)+1,j),j=1,3),(node_list(3*(i-1)+2,j),j=1,3)
2421 
2422 
2423  ENDIF
2424  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 2) THEN
2425 
2426  READ(11,*) (node_list(3*(i-1)+1,j),j=1,3)
2427 
2428  ENDIF
2429 
2430 
2431 ! EMMNLEMNTTMENLMNTMT load data for ELEMENT CONNECTIVITY info
2432 
2433  READ(11,*) string
2434 
2435  DO i=1,total_number_of_elements
2436 
2437  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
2438 
2439 ! we have nodes starting number of 0
2440  DO j=1,number_of_nodes_per_element
2441  element_list(i,j)=element_list(i,j)+1
2442  ENDDO
2443 
2444  DO j=1,number_of_nodes_per_element-1
2445 
2446  DO k=1+j,number_of_nodes_per_element
2447  there_is_in_connectivity_list = 0
2448  DO n=1,connectivity_number(element_list(i,j))
2449 ! print *,"elem",I,"connectivity number",N
2450 
2451  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
2452  there_is_in_connectivity_list = 1
2453  ENDIF
2454 
2455  ENDDO
2456 
2457  IF (there_is_in_connectivity_list .EQ. 0) THEN
2458 
2459  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
2460  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
2461  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
2462  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
2463 
2464  ENDIF
2465 
2466  ENDDO
2467 
2468  ENDDO
2469 
2470  ENDDO
2471 
2472  CLOSE(11)
2473 
2474 ! SDEESDSEDSEESDSEEDS load SEED VALUES data at the nods
2475 ! set input for *seed points* = LIST
2476  IF (input_type_for_seed_value .EQ. "LIST") THEN
2477 
2478  DO i=1,total_number_of_nodes
2479 
2480  IF (status_mask(i) .NE. "SEED POINT") THEN
2481  seed_value(i) = 1000.0_dp
2482  ENDIF
2483 
2484  ENDDO
2485 
2486  ENDIF
2487 
2488 ! read input for *seed points* = FILE
2489  IF (input_type_for_seed_value .EQ. "FILE") THEN
2490 
2491  text_length = index(input_file_name,' ') - 1
2492  string = input_file_name(1:text_length)//".estm"
2493 
2494  DO i=1,total_number_of_nodes
2495  status_mask(i) = ""
2496  ENDDO
2497 
2498  OPEN (11,file=string)
2499  READ (11,*) n
2500 
2501  DO i=1,n
2502 
2503  READ(11,*) j,seed_value(j+1)
2504 ! PRINT*,(NODE_LIST(J+1,K),K=1,3),SEED_VALUE(J+1)
2505  status_mask(j+1) = "SEED POINT"
2506 
2507  ENDDO
2508 
2509  CLOSE(11)
2510 
2511  DO i=1,total_number_of_nodes
2512 
2513  IF (status_mask(i) .NE. "SEED POINT") THEN
2514  seed_value(i) = 1000.0_dp
2515  ENDIF
2516 
2517  ENDDO
2518 
2519  ENDIF
2520 
2521 ! NDCONDNCODNCODNDCOM load CONDUCTIVITY TENSOR data for the nods
2522 ! for the ISOTROPIC materials
2523  IF (material_behaviour .EQ. "ISOTROPIC") THEN
2524 
2525  DO i=1,total_number_of_nodes
2526  conductivity_tensor(i,1)=1.0_dp
2527  conductivity_tensor(i,2)=0.0_dp
2528  conductivity_tensor(i,3)=0.0_dp
2529  conductivity_tensor(i,4)=0.0_dp
2530  conductivity_tensor(i,5)=1.0_dp
2531  conductivity_tensor(i,6)=0.0_dp
2532  conductivity_tensor(i,7)=0.0_dp
2533  conductivity_tensor(i,8)=0.0_dp
2534  conductivity_tensor(i,9)=1.0_dp
2535  ENDDO
2536 
2537  ENDIF
2538 
2539 ! for the ANISOTROPIC materials
2540  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
2541 
2542  text_length = index(input_file_name,' ') - 1
2543  string = input_file_name(1:text_length)//".fiber"
2544 
2545  OPEN (11,file=string)
2546  READ (11,*) string
2547 
2548  IF (input_type_for_conductivity .EQ. "TENSOR") THEN
2549  DO i=1,total_number_of_nodes
2550  READ(11,*) string,(conductivity_tensor(i,j),j=1,9)
2551  ENDDO
2552  ENDIF
2553 
2554  IF (input_type_for_conductivity .EQ. "VECTOR") THEN
2555  IF (string .EQ. '#') then! begin if
2556 
2557  DO WHILE (string .NE. 'fiber')
2558 
2559  READ(11,*) string
2560 
2561  ENDDO
2562 
2563 ! PRINT *,STRING
2564 
2565  DO i=1,int((total_number_of_nodes/3.0_dp)+0.5_dp)-1
2566 
2567  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3),(conductivity_tensor(3*(i-1)+2,j),j=1,3), &
2568  & (conductivity_tensor(3*(i-1)+3,j),j=1,3)
2569 
2570  ENDDO
2571 
2572  i=int((total_number_of_nodes/3.0_dp)+0.5_dp)
2573 
2574  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 0) THEN
2575 
2576  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3),(conductivity_tensor(3*(i-1)+2,j),j=1,3), &
2577  & (conductivity_tensor(3*(i-1)+3,j),j=1,3)
2578 
2579  ENDIF
2580  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 1) THEN
2581 
2582  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3),(conductivity_tensor(3*(i-1)+2,j),j=1,3)
2583 
2584  ENDIF
2585  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 2) THEN
2586 
2587  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3)
2588 
2589  ENDIF
2590 
2591 
2592  ELSE
2593  DO i=1,total_number_of_nodes
2594 
2595  READ(11,*) (conductivity_tensor(i,j),j=1,3)
2596 
2597  ENDDO
2598  ENDIF ! end if
2599 
2600  DO i=1,total_number_of_nodes
2601 
2602 ! READ(11,*) (CONDUCTIVITY_TENSOR(I,J),J=1,3)
2603 
2604  a=(/conductivity_tensor(i,1),conductivity_tensor(i,2),conductivity_tensor(i,3)/)
2605 
2606  b=(/0.0_dp,0.0_dp,1.0_dp/)
2607  CALL vector_vector_product(a,b,dot_product_value,err)
2608 
2609  IF (dot_product_value .LT. 0.8_dp) THEN
2610  CALL cross_product(a,b,c,err,error,*999)
2611  ELSE
2612  b=(/0.0_dp,1.0_dp,0.0_dp/)
2613  CALL vector_vector_product(a,b,dot_product_value,err)
2614  IF (dot_product_value .LT. 0.8_dp) THEN
2615  CALL cross_product(a,b,c,err,error,*999)
2616  ELSE
2617  b=(/1.0_dp,0.0_dp,0.0_dp/)
2618  CALL cross_product(a,b,c,err,error,*999)
2619  ENDIF
2620  ENDIF
2621 
2622  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>zero_tolerance) THEN
2623 
2624  conductivity_tensor(i,4) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2625  conductivity_tensor(i,5) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2626  conductivity_tensor(i,6) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2627 
2628  ENDIF
2629 
2630  b=(/conductivity_tensor(i,4),conductivity_tensor(i,5),conductivity_tensor(i,6)/)
2631  CALL cross_product(a,b,c,err,error,*999)
2632 
2633  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>=zero_tolerance) THEN
2634 
2635  conductivity_tensor(i,7) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2636  conductivity_tensor(i,8) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2637  conductivity_tensor(i,9) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2638 
2639  ENDIF
2640 
2641 ! PRINT*,(CONDUCTIVITY_TENSOR(I,J),J=1,3),DOT_PRODUCT_VALUE
2642 
2643  ENDDO
2644  ENDIF
2645 
2646  CLOSE(11)
2647 
2648  ENDIF
2649 
2650 ! FNCONSJDCNFUCNSUCNF load VELOCITY FUNCTION data for the nods
2651 ! for the ISOTROPIC materials
2652  IF (material_behaviour .EQ. "ISOTROPIC") THEN
2653 
2654 ! if it comes with FIXED format
2655  IF (input_type_for_speed_function .EQ. "FIXED") THEN
2656 
2657  DO i=1,total_number_of_nodes
2658  DO j=1,3
2659  speed_function_table(i,j) = speed_function_along_eigen_vector(1)
2660  ENDDO
2661  ENDDO
2662 
2663  ENDIF
2664 
2665 ! if it comes with FILE format
2666  IF (input_type_for_speed_function .EQ. "FILE") THEN
2667 
2668  text_length = index(input_file_name,' ') - 1
2669  string = input_file_name(1:text_length)//".cond"
2670 
2671  OPEN (11,file=string)
2672  READ (11,*) n
2673 
2674  DO i=1,n
2675 
2676  READ(11,*) j,speed_function_table(j,1)
2677  speed_function_table(j,2) = speed_function_table(j,1)
2678  speed_function_table(j,3) = speed_function_table(j,1)
2679 
2680  ENDDO
2681 
2682  CLOSE(11)
2683 
2684  ENDIF
2685 
2686  ENDIF
2687 
2688 ! for the ANISOTROPIC materials
2689  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
2690 
2691 ! if it comes with FIXED format
2692  IF (input_type_for_speed_function .EQ. "FIXED") THEN
2693 
2694  DO i=1,total_number_of_nodes
2695  DO j=1,3
2696  speed_function_table(i,j) = 0.0_dp
2697  ENDDO
2698  ENDDO
2699 
2700  DO i=1,total_number_of_nodes
2701 ! DO J=1,CONNECTIVITY_NUMBER(I)
2702 ! IF (CONNECTIVITY_LIST(I,J) > TOTAL_NUMBER_OF_NODES) THEN
2703 ! PRINT*, I,J,CONNECTIVITY_NUMBER(I),CONNECTIVITY_LIST(I,J),TOTAL_NUMBER_OF_NODES
2704 ! ENDIF
2705  speed_function_table(i,1) = speed_function_along_eigen_vector(1)
2706  speed_function_table(i,2) = speed_function_along_eigen_vector(2)
2707  speed_function_table(i,3) = speed_function_along_eigen_vector(3)
2708 ! SPEED_FUNCTION_TABLE(I,CONNECTIVITY_LIST(I,J),1) = SPEED_FUNCTION_ALONG_EIGEN_VECTOR(1)
2709 ! SPEED_FUNCTION_TABLE(I,CONNECTIVITY_LIST(I,J),2) = SPEED_FUNCTION_ALONG_EIGEN_VECTOR(2)
2710 ! SPEED_FUNCTION_TABLE(I,CONNECTIVITY_LIST(I,J),3) = SPEED_FUNCTION_ALONG_EIGEN_VECTOR(3)
2711 ! ENDDO
2712  ENDDO
2713 
2714  ENDIF
2715 
2716 ! if it comes with FILE format
2717  IF (input_type_for_speed_function .EQ. "FILE") THEN
2718 
2719  text_length = index(input_file_name,' ') - 1
2720 
2721  string = input_file_name(1:text_length)//".cond"
2722 
2723  OPEN (11,file=string)
2724  READ (11,*) n
2725 
2726  DO i=1,n
2727 
2728  READ(11,*) j,speed_function_table(j,1),speed_function_table(j,2),speed_function_table(j,3)
2729 
2730  ENDDO
2731 
2732  CLOSE(11)
2733 
2734  ENDIF
2735 
2736  ENDIF
2737 
2738  ENDIF
2739 ! """""""""""""""""""""""""""""""""""INPUT OF VTK FORMAT""""""""""""""""""""""""""""""
2740  IF (input_file_format .EQ. "VTKTET1NPL") THEN
2741 
2742  number_of_nodes_per_element = 4
2743 
2744 ! NDNAPDNDONOEEENODED load data for NODAL POSITIONING info
2745  text_length = index(input_file_name,' ') - 1
2746  string = input_file_name(1:text_length)//".vtk"
2747 
2748  OPEN (11,file=string)
2749  READ(11,*) string
2750  READ(11,*) string
2751  READ(11,*) string
2752  READ(11,*) string
2753  READ(11,*) string
2754 
2755  DO i=1,total_number_of_nodes
2756 
2757  READ(11,*) (node_list(i,j),j=1,3)
2758 
2759  ENDDO
2760 
2761 ! EMMNLEMNTTMENLMNTMT load data for ELEMENT CONNECTIVITY info
2762 
2763  READ(11,*) string
2764 
2765  DO n=1,total_number_of_nodes
2766  connectivity_number(n)=0
2767  ENDDO
2768 
2769  DO i=1,total_number_of_elements
2770 
2771  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
2772 
2773 ! we have nodes starting number of 0
2774  DO j=1,number_of_nodes_per_element
2775  element_list(i,j)=element_list(i,j)+1
2776  ENDDO
2777 
2778  DO j=1,number_of_nodes_per_element-1
2779 
2780  DO k=1+j,number_of_nodes_per_element
2781  there_is_in_connectivity_list = 0
2782  DO n=1,connectivity_number(element_list(i,j))
2783 ! print *,"elem",I,"connectivity number",N
2784 
2785  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
2786  there_is_in_connectivity_list = 1
2787  ENDIF
2788 
2789  ENDDO
2790 
2791  IF (there_is_in_connectivity_list .EQ. 0) THEN
2792 
2793  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
2794  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
2795  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
2796  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
2797 
2798  ENDIF
2799 
2800  ENDDO
2801 
2802  ENDDO
2803 
2804  ENDDO
2805 
2806  CLOSE(11)
2807 
2808 ! SDEESDSEDSEESDSEEDS load SEED VALUES data at the nods
2809 ! set input for *seed points* = LIST
2810  IF (input_type_for_seed_value .EQ. "LIST") THEN
2811 
2812  DO i=1,total_number_of_nodes
2813 
2814  IF (status_mask(i) .NE. "SEED POINT") THEN
2815  seed_value(i) = 1000.0_dp
2816  ENDIF
2817 
2818  ENDDO
2819 
2820  ENDIF
2821 
2822 ! read input for *seed points* = FILE
2823  IF (input_type_for_seed_value .EQ. "FILE") THEN
2824 
2825  text_length = index(input_file_name,' ') - 1
2826  string = input_file_name(1:text_length)//".estm"
2827 
2828  OPEN (11,file=string)
2829  READ (11,*) n
2830 
2831  DO i=1,n
2832 
2833  READ(11,*) j,seed_value(j+1)
2834 
2835  status_mask(j+1) = "SEED POINT"
2836 
2837  ENDDO
2838 
2839  CLOSE(11)
2840 
2841  DO i=1,total_number_of_nodes
2842 
2843  IF (status_mask(i) .NE. "SEED POINT") THEN
2844  seed_value(i) = 1000.0_dp
2845  ENDIF
2846 
2847  ENDDO
2848 
2849  ENDIF
2850 
2851 ! NDCONDNCODNCODNDCOM load CONDUCTIVITY TENSOR data for the nods
2852 ! for the ISOTROPIC materials
2853  IF (material_behaviour .EQ. "ISOTROPIC") THEN
2854 
2855  DO i=1,total_number_of_nodes
2856  conductivity_tensor(i,1)=1.0_dp
2857  conductivity_tensor(i,2)=0.0_dp
2858  conductivity_tensor(i,3)=0.0_dp
2859  conductivity_tensor(i,4)=0.0_dp
2860  conductivity_tensor(i,5)=1.0_dp
2861  conductivity_tensor(i,6)=0.0_dp
2862  conductivity_tensor(i,7)=0.0_dp
2863  conductivity_tensor(i,8)=0.0_dp
2864  conductivity_tensor(i,9)=1.0_dp
2865  ENDDO
2866 
2867  ENDIF
2868 
2869 ! for the ANISOTROPIC materials
2870  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
2871 
2872  text_length = index(input_file_name,' ') - 1
2873  string = input_file_name(1:text_length)//".fiber"
2874 
2875  OPEN (11,file=string)
2876  READ (11,*) string
2877 ! print *,STRING
2878 
2879  IF (input_type_for_conductivity .EQ. "TENSOR") THEN
2880  DO i=1,total_number_of_nodes
2881  READ(11,*) string,(conductivity_tensor(i,j),j=1,9)
2882  ENDDO
2883  ENDIF
2884 
2885  IF (input_type_for_conductivity .EQ. "VECTOR") THEN
2886  IF (string .EQ. '#') then! begin if
2887 
2888  DO WHILE (string .NE. 'fiber')
2889 
2890  READ(11,*) string
2891 
2892  ENDDO
2893 
2894 ! PRINT *,STRING
2895 
2896  DO i=1,int((total_number_of_nodes/3.0_dp)+0.5_dp)-1
2897 
2898  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3),(conductivity_tensor(3*(i-1)+2,j),j=1,3), &
2899  & (conductivity_tensor(3*(i-1)+3,j),j=1,3)
2900 
2901  ENDDO
2902 
2903  i=int((total_number_of_nodes/3.0_dp)+0.5_dp)
2904 
2905  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 0) THEN
2906 
2907  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3),(conductivity_tensor(3*(i-1)+2,j),j=1,3), &
2908  & (conductivity_tensor(3*(i-1)+3,j),j=1,3)
2909 
2910  ENDIF
2911  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 1) THEN
2912 
2913  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3),(conductivity_tensor(3*(i-1)+2,j),j=1,3)
2914 
2915  ENDIF
2916  IF (3*int((total_number_of_nodes/3.0_dp)+0.5_dp)-total_number_of_nodes .EQ. 2) THEN
2917 
2918  READ(11,*) (conductivity_tensor(3*(i-1)+1,j),j=1,3)
2919 
2920  ENDIF
2921 
2922 
2923  ELSE
2924  DO i=1,total_number_of_nodes
2925 
2926  READ(11,*) (conductivity_tensor(i,j),j=1,3)
2927 
2928  ENDDO
2929  ENDIF ! end if
2930 
2931  DO i=1,total_number_of_nodes
2932 
2933 ! READ(11,*) (CONDUCTIVITY_TENSOR(I,J),J=1,3)
2934 
2935  a=(/conductivity_tensor(i,1),conductivity_tensor(i,2),conductivity_tensor(i,3)/)
2936  b=(/0.0_dp,0.0_dp,1.0_dp/)
2937  CALL vector_vector_product(a,b,dot_product_value,err)
2938 
2939  IF (dot_product_value .LT. 0.8_dp) THEN
2940  CALL cross_product(a,b,c,err,error,*999)
2941  ELSE
2942  b=(/0.0_dp,1.0_dp,0.0_dp/)
2943  CALL vector_vector_product(a,b,dot_product_value,err)
2944  IF (dot_product_value .LT. 0.8_dp) THEN
2945  CALL cross_product(a,b,c,err,error,*999)
2946  ELSE
2947  b=(/1.0_dp,0.0_dp,0.0_dp/)
2948  CALL cross_product(a,b,c,err,error,*999)
2949  ENDIF
2950  ENDIF
2951 
2952  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>=zero_tolerance) THEN
2953 
2954  conductivity_tensor(i,4) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2955  conductivity_tensor(i,5) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2956  conductivity_tensor(i,6) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2957 
2958  ENDIF
2959 
2960  b=(/conductivity_tensor(i,4),conductivity_tensor(i,5),conductivity_tensor(i,6)/)
2961  CALL cross_product(a,b,c,err,error,*999)
2962 
2963  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>=zero_tolerance) THEN
2964 
2965  conductivity_tensor(i,7) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2966  conductivity_tensor(i,8) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2967  conductivity_tensor(i,9) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
2968 
2969  ENDIF
2970 
2971 ! PRINT*,(CONDUCTIVITY_TENSOR(I,J),J=1,3),DOT_PRODUCT_VALUE
2972 
2973  ENDDO
2974  ENDIF
2975 
2976  CLOSE(11)
2977 
2978  ENDIF
2979 
2980 ! FNCONSJDCNFUCNSUCNF load VELOCITY FUNCTION data for the nods
2981 ! for the ISOTROPIC materials
2982  IF (material_behaviour .EQ. "ISOTROPIC") THEN
2983 
2984 ! if it comes with FIXED format
2985  IF (input_type_for_speed_function .EQ. "FIXED") THEN
2986 
2987  DO i=1,total_number_of_nodes
2988  DO j=1,3
2989  speed_function_table(i,j) = speed_function_along_eigen_vector(1)
2990  ENDDO
2991  ENDDO
2992 
2993  ENDIF
2994 
2995 ! if it comes with FILE format
2996  IF (input_type_for_speed_function .EQ. "FILE") THEN
2997 
2998  text_length = index(input_file_name,' ') - 1
2999  string = input_file_name(1:text_length)//".cond"
3000 
3001  OPEN (11,file=string)
3002  READ (11,*) n
3003 
3004  DO i=1,n
3005 
3006  READ(11,*) j,speed_function_table(j,1)
3007  speed_function_table(j,2) = speed_function_table(j,1)
3008  speed_function_table(j,3) = speed_function_table(j,1)
3009 
3010  ENDDO
3011 
3012  CLOSE(11)
3013 
3014  ENDIF
3015 
3016  ENDIF
3017 
3018 ! for the ANISOTROPIC materials
3019  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
3020 
3021 ! if it comes with FIXED format
3022  IF (input_type_for_speed_function .EQ. "FIXED") THEN
3023 
3024  DO i=1,total_number_of_nodes
3025  DO j=1,3
3026  speed_function_table(i,j) = speed_function_along_eigen_vector(j)
3027  ENDDO
3028  ENDDO
3029 
3030  ENDIF
3031 
3032 ! if it comes with FILE format
3033  IF (input_type_for_speed_function .EQ. "FILE") THEN
3034 
3035  text_length = index(input_file_name,' ') - 1
3036 
3037  string = input_file_name(1:text_length)//".cond"
3038 
3039  OPEN (11,file=string)
3040  READ (11,*) n
3041 
3042  DO i=1,n
3043 
3044  READ(11,*) j,speed_function_table(j,1),speed_function_table(j,2),&
3045  &speed_function_table(j,3)
3046 
3047  ENDDO
3048 
3049  CLOSE(11)
3050 
3051  ENDIF
3052 
3053  ENDIF
3054 
3055  ENDIF
3056 
3057 ! """""""""""""""""""""""""""""""""""INPUT OF CARP FORMAT""""""""""""""""""""""""""""""
3058  IF (input_file_format .EQ. "CARP") THEN
3059 
3060  number_of_nodes_per_element = 4
3061 
3062 ! NDNAPDNDONOEEENODED load data for NODAL POSITIONING info
3063  text_length = index(input_file_name,' ') - 1
3064  string = input_file_name(1:text_length)//".pts"
3065 
3066  OPEN (11,file=string)
3067  READ (11,*) string
3068 
3069  DO i=1,total_number_of_nodes
3070 
3071  READ(11,*) (node_list(i,j),j=1,3)
3072 
3073  ENDDO
3074 
3075  CLOSE(11)
3076 
3077 ! EMMNLEMNTTMENLMNTMT load data for ELEMENT CONNECTIVITY info
3078  text_length = index(input_file_name,' ') - 1
3079  string = input_file_name(1:text_length)//".elem"
3080 
3081  OPEN (11,file=string)
3082  READ (11,*) total_number_of_elements
3083 
3084  DO n=1,total_number_of_nodes
3085  connectivity_number(n)=0
3086  ENDDO
3087 
3088  DO i=1,total_number_of_elements
3089 
3090  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
3091 
3092  DO j=1,number_of_nodes_per_element
3093  element_list(i,j)=element_list(i,j)+1
3094  ENDDO
3095 
3096  DO j=1,number_of_nodes_per_element-1
3097 
3098  DO k=1+j,number_of_nodes_per_element
3099  there_is_in_connectivity_list = 0
3100  DO n=1,connectivity_number(element_list(i,j))
3101 ! print *,"elem",I,"connectivity number",N
3102 
3103  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
3104  there_is_in_connectivity_list = 1
3105  ENDIF
3106 
3107  ENDDO
3108 
3109  IF (there_is_in_connectivity_list .EQ. 0) THEN
3110 
3111  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
3112  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
3113  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
3114  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
3115 
3116  ENDIF
3117 
3118  ENDDO
3119 
3120  ENDDO
3121 
3122  ENDDO
3123 
3124  CLOSE(11)
3125 
3126 ! SDEESDSEDSEESDSEEDS load SEED VALUES data at the nods
3127 ! set input for *seed points* = LIST
3128  IF (input_type_for_seed_value .EQ. "LIST") THEN
3129 
3130  DO i=1,total_number_of_nodes
3131 
3132  IF (status_mask(i) .NE. "SEED POINT") THEN
3133  seed_value(i) = 1000.0_dp
3134  ENDIF
3135 
3136  ENDDO
3137 
3138  ENDIF
3139 
3140 ! read input for *seed points* = FILE
3141  IF (input_type_for_seed_value .EQ. "FILE") THEN
3142 
3143  text_length = index(input_file_name,' ') - 1
3144  string = input_file_name(1:text_length)//".estm"
3145 
3146  OPEN (11,file=string)
3147  READ (11,*) string
3148 
3149  DO i=1,total_number_of_nodes
3150 
3151  READ(11,*) string,seed_value(i)
3152 
3153  ENDDO
3154 
3155  CLOSE(11)
3156 
3157  ENDIF
3158 
3159 ! NDCONDNCODNCODNDCOM load CONDUCTIVITY TENSOR data for the nods
3160 ! for the ISOTROPIC materials
3161  IF (material_behaviour .EQ. "ISOTROPIC") THEN
3162 
3163  DO i=1,total_number_of_nodes
3164  conductivity_tensor(i,1)=1.0_dp
3165  conductivity_tensor(i,2)=0.0_dp
3166  conductivity_tensor(i,3)=0.0_dp
3167  conductivity_tensor(i,4)=0.0_dp
3168  conductivity_tensor(i,5)=1.0_dp
3169  conductivity_tensor(i,6)=0.0_dp
3170  conductivity_tensor(i,7)=0.0_dp
3171  conductivity_tensor(i,8)=0.0_dp
3172  conductivity_tensor(i,9)=1.0_dp
3173  ENDDO
3174 
3175  ENDIF
3176 
3177 ! for the ANISOTROPIC materials
3178  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
3179 
3180  text_length = index(input_file_name,' ') - 1
3181  string = input_file_name(1:text_length)//".lon"
3182 
3183  OPEN (11,file=string)
3184  READ (11,*) string
3185 
3186  IF (input_type_for_conductivity .EQ. "TENSOR") THEN
3187  DO i=1,total_number_of_nodes
3188  READ(11,*) string,(conductivity_tensor(i,j),j=1,9)
3189  ENDDO
3190  ENDIF
3191 
3192  IF (input_type_for_conductivity .EQ. "VECTOR" .OR. input_type_for_conductivity .EQ. " ") THEN
3193  DO i=1,total_number_of_elements
3194  READ(11,*) (conductivity_tensor(element_list(i,1),j),j=1,3)
3195 
3196  a=(/conductivity_tensor(element_list(i,1),1),conductivity_tensor(element_list(i,1),2), &
3197  & conductivity_tensor(element_list(i,1),3)/)
3198  b=(/0.0_dp,0.0_dp,1.0_dp/)
3199  CALL vector_vector_product(a,b,dot_product_value,err)
3200 
3201  IF (dot_product_value .LT. 0.8_dp) THEN
3202  CALL cross_product(a,b,c,err,error,*999)
3203  ELSE
3204  b=(/0.0_dp,1.0_dp,0.0_dp/)
3205  CALL vector_vector_product(a,b,dot_product_value,err)
3206  IF (dot_product_value .LT. 0.8_dp) THEN
3207  CALL cross_product(a,b,c,err,error,*999)
3208  ELSE
3209  b=(/1.0_dp,0.0_dp,0.0_dp/)
3210  CALL cross_product(a,b,c,err,error,*999)
3211  ENDIF
3212  ENDIF
3213 
3214  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>=zero_tolerance) THEN
3215 
3216  conductivity_tensor(element_list(i,1),4) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3217  conductivity_tensor(element_list(i,1),5) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3218  conductivity_tensor(element_list(i,1),6) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3219 
3220  ENDIF
3221 
3222  b=(/conductivity_tensor(element_list(i,1),4),conductivity_tensor(element_list(i,1),5), &
3223  & conductivity_tensor(element_list(i,1),6)/)
3224  CALL cross_product(a,b,c,err,error,*999)
3225 
3226  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>=zero_tolerance) THEN
3227 
3228  conductivity_tensor(element_list(i,1),7) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3229  conductivity_tensor(element_list(i,1),8) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3230  conductivity_tensor(element_list(i,1),9) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3231 
3232  ENDIF
3233 
3234  DO j=2,number_of_nodes_per_element
3235  DO k=1,9
3236  conductivity_tensor(element_list(i,j),k)=conductivity_tensor(element_list(i,1),k)
3237  ENDDO
3238  ENDDO
3239 
3240 ! PRINT*,(CONDUCTIVITY_TENSOR(I,J),J=1,3),DOT_PRODUCT_VALUE
3241 
3242  ENDDO
3243  ENDIF
3244 
3245  CLOSE(11)
3246 
3247  ENDIF
3248 
3249 ! FNCONSJDCNFUCNSUCNF load VELOCITY FUNCTION data for the nods
3250 ! for the ISOTROPIC materials
3251  IF (material_behaviour .EQ. "ISOTROPIC") THEN
3252 
3253 ! if it comes with FIXED format
3254  IF (input_type_for_speed_function .EQ. "FIXED") THEN
3255 
3256  DO i=1,total_number_of_nodes
3257  DO j=1,3
3258  speed_function_table(i,j) = speed_function_along_eigen_vector(1)
3259  ENDDO
3260  ENDDO
3261 
3262  ENDIF
3263 
3264 ! if it comes with FILE format
3265  IF (input_type_for_speed_function .EQ. "FILE") THEN
3266 
3267  text_length = index(input_file_name,' ') - 1
3268  string = input_file_name(1:text_length)//".cond"
3269 
3270  OPEN (11,file=string)
3271  READ (11,*) n
3272 
3273  DO i=1,n
3274 
3275  READ(11,*) j,speed_function_table(j,1)
3276  speed_function_table(j,2) = speed_function_table(j,1)
3277  speed_function_table(j,3) = speed_function_table(j,1)
3278 
3279  ENDDO
3280 
3281  CLOSE(11)
3282 
3283  ENDIF
3284 
3285  ENDIF
3286 
3287 ! for the ANISOTROPIC materials
3288  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
3289 
3290 ! if it comes with FIXED format
3291  IF (input_type_for_speed_function .EQ. "FIXED") THEN
3292 
3293  DO i=1,total_number_of_nodes
3294  DO j=1,3
3295  speed_function_table(i,j) = speed_function_along_eigen_vector(j)
3296  ENDDO
3297  ENDDO
3298 
3299  ENDIF
3300 
3301 ! if it comes with FILE format
3302  IF (input_type_for_speed_function .EQ. "FILE") THEN
3303 
3304  text_length = index(input_file_name,' ') - 1
3305 
3306  string = input_file_name(1:text_length)//".cond"
3307 
3308  OPEN (11,file=string)
3309  READ (11,*) n
3310 
3311  DO i=1,n
3312 
3313  READ(11,*) j,speed_function_table(j,1),speed_function_table(j,2),&
3314  &speed_function_table(j,3)
3315 
3316  ENDDO
3317 
3318  CLOSE(11)
3319 
3320  ENDIF
3321 
3322  ENDIF
3323 
3324  ENDIF
3325 
3326 ! """""""""""""""""""""""""""""""""""INPUT OF TETGEN FORMAT""""""""""""""""""""""""""""""
3327  IF (input_file_format .EQ. "TETGEN") THEN
3328 
3329  number_of_nodes_per_element = 4
3330 
3331 ! NDNAPDNDONOEEENODED load data for NODAL POSITIONING info
3332  text_length = index(input_file_name,' ') - 1
3333  string = input_file_name(1:text_length)//".node"
3334 
3335  OPEN (11,file=string)
3336  READ (11,*) string
3337 
3338  READ(11,*) first_node_number,(node_list(1,j),j=1,3)
3339 
3340  DO i=2,total_number_of_nodes
3341 
3342 
3343  READ(11,*) string,(node_list(i,j),j=1,3)
3344 
3345 
3346  ENDDO
3347 
3348  CLOSE(11)
3349 
3350 ! EMMNLEMNTTMENLMNTMT load data for ELEMENT CONNECTIVITY info
3351  text_length = index(input_file_name,' ') - 1
3352  string = input_file_name(1:text_length)//".ele"
3353 
3354  OPEN (11,file=string)
3355  READ (11,*) total_number_of_elements
3356 
3357  DO n=1,total_number_of_nodes
3358  connectivity_number(n)=0
3359  ENDDO
3360 
3361  DO i=1,total_number_of_elements
3362 
3363  READ(11,*) string,(element_list(i,j),j=1,number_of_nodes_per_element)
3364 
3365  IF (first_node_number .EQ. 0) THEN
3366  DO j=1,number_of_nodes_per_element
3367  element_list(i,j)=element_list(i,j)+1
3368  ENDDO
3369  ENDIF
3370 
3371  DO j=1,number_of_nodes_per_element-1
3372 
3373  DO k=1+j,number_of_nodes_per_element
3374  there_is_in_connectivity_list = 0
3375  DO n=1,connectivity_number(element_list(i,j))
3376 ! print *,"elem",I,"connectivity number",N
3377 
3378  IF (connectivity_list(element_list(i,j),n) .EQ. element_list(i,k)) THEN
3379  there_is_in_connectivity_list = 1
3380  ENDIF
3381 
3382  ENDDO
3383 
3384  IF (there_is_in_connectivity_list .EQ. 0) THEN
3385 
3386  connectivity_number(element_list(i,j)) = connectivity_number(element_list(i,j)) + 1
3387  connectivity_list(element_list(i,j),connectivity_number(element_list(i,j))) = element_list(i,k)
3388  connectivity_number(element_list(i,k)) = connectivity_number(element_list(i,k)) + 1
3389  connectivity_list(element_list(i,k),connectivity_number(element_list(i,k))) = element_list(i,j)
3390 
3391  ENDIF
3392 
3393  ENDDO
3394 
3395  ENDDO
3396 
3397  ENDDO
3398 
3399  CLOSE(11)
3400 
3401 ! SDEESDSEDSEESDSEEDS load SEED VALUES data at the nods
3402 ! set input for *seed points* = LIST
3403  IF (input_type_for_seed_value .EQ. "LIST") THEN
3404 
3405  DO i=1,total_number_of_nodes
3406 
3407  IF (status_mask(i) .NE. "SEED POINT") THEN
3408  seed_value(i) = 1000.0_dp
3409  ENDIF
3410 
3411  ENDDO
3412 
3413  ENDIF
3414 
3415 ! read input for *seed points* = FILE
3416  IF (input_type_for_seed_value .EQ. "FILE") THEN
3417 
3418  text_length = index(input_file_name,' ') - 1
3419  string = input_file_name(1:text_length)//".estm"
3420 
3421  OPEN (11,file=string)
3422  READ (11,*) string
3423 
3424  DO i=1,total_number_of_nodes
3425 
3426  READ(11,*) string,seed_value(i)
3427 
3428  ENDDO
3429 
3430  CLOSE(11)
3431 
3432  ENDIF
3433 
3434 ! NDCONDNCODNCODNDCOM load CONDUCTIVITY TENSOR data for the nods
3435 ! for the ISOTROPIC materials
3436  IF (material_behaviour .EQ. "ISOTROPIC") THEN
3437 
3438  DO i=1,total_number_of_nodes
3439  conductivity_tensor(i,1)=1.0_dp
3440  conductivity_tensor(i,2)=0.0_dp
3441  conductivity_tensor(i,3)=0.0_dp
3442  conductivity_tensor(i,4)=0.0_dp
3443  conductivity_tensor(i,5)=1.0_dp
3444  conductivity_tensor(i,6)=0.0_dp
3445  conductivity_tensor(i,7)=0.0_dp
3446  conductivity_tensor(i,8)=0.0_dp
3447  conductivity_tensor(i,9)=1.0_dp
3448  ENDDO
3449 
3450  ENDIF
3451 
3452 ! for the ANISOTROPIC materials
3453  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
3454 
3455  text_length = index(input_file_name,' ') - 1
3456  string = input_file_name(1:text_length)//".fiber"
3457 
3458  OPEN (11,file=string)
3459  READ (11,*) string
3460 
3461  IF (input_type_for_conductivity .EQ. "TENSOR") THEN
3462  DO i=1,total_number_of_nodes
3463  READ(11,*) string,(conductivity_tensor(i,j),j=1,9)
3464  ENDDO
3465  ENDIF
3466 
3467  IF (input_type_for_conductivity .EQ. "VECTOR") THEN
3468  DO i=1,total_number_of_nodes
3469  READ(11,*) string,(conductivity_tensor(i,j),j=1,3)
3470 
3471  a=(/conductivity_tensor(i,1),conductivity_tensor(i,2),conductivity_tensor(i,3)/)
3472  b=(/0.0_dp,0.0_dp,1.0_dp/)
3473  CALL vector_vector_product(a,b,dot_product_value,err)
3474 
3475  IF (dot_product_value .LT. 0.8_dp) THEN
3476  CALL cross_product(a,b,c,err,error,*999)
3477  ELSE
3478  b=(/0.0_dp,1.0_dp,0.0_dp/)
3479  CALL vector_vector_product(a,b,dot_product_value,err)
3480  IF (dot_product_value .LT. 0.8_dp) THEN
3481  CALL cross_product(a,b,c,err,error,*999)
3482  ELSE
3483 
3484  b=(/1.0_dp,0.0_dp,0.0_dp/)
3485  CALL cross_product(a,b,c,err,error,*999)
3486  ENDIF
3487  ENDIF
3488 
3489  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>zero_tolerance) THEN
3490 
3491  conductivity_tensor(i,4) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3492  conductivity_tensor(i,5) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3493  conductivity_tensor(i,6) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3494 
3495  ENDIF
3496 
3497  b=(/conductivity_tensor(i,4),conductivity_tensor(i,5),conductivity_tensor(i,6)/)
3498  CALL cross_product(a,b,c,err,error,*999)
3499 
3500  IF (abs(sqrt(c(1)**2+c(2)**2+c(3)**2))>zero_tolerance) THEN
3501 
3502  conductivity_tensor(i,7) = c(1)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3503  conductivity_tensor(i,8) = c(2)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3504  conductivity_tensor(i,9) = c(3)/sqrt(c(1)**2+c(2)**2+c(3)**2)
3505 
3506  ENDIF
3507 
3508 ! PRINT*,(CONDUCTIVITY_TENSOR(I,J),J=1,3),DOT_PRODUCT_VALUE
3509 
3510  ENDDO
3511  ENDIF
3512 
3513  CLOSE(11)
3514 
3515  ENDIF
3516 
3517 ! FNCONSJDCNFUCNSUCNF load VELOCITY FUNCTION data for the nods
3518 ! for the ISOTROPIC materials
3519  IF (material_behaviour .EQ. "ISOTROPIC") THEN
3520 
3521 ! if it comes with FIXED format
3522  IF (input_type_for_speed_function .EQ. "FIXED") THEN
3523 
3524  DO i=1,total_number_of_nodes
3525  DO j=1,3
3526  speed_function_table(i,j) = speed_function_along_eigen_vector(1)
3527  ENDDO
3528  ENDDO
3529 
3530  ENDIF
3531 
3532 ! if it comes with FILE format
3533  IF (input_type_for_speed_function .EQ. "FILE") THEN
3534 
3535  text_length = index(input_file_name,' ') - 1
3536  string = input_file_name(1:text_length)//".cond"
3537 
3538  OPEN (11,file=string)
3539  READ (11,*) n
3540 
3541  DO i=1,n
3542 
3543  READ(11,*) j,speed_function_table(j,1)
3544  speed_function_table(j,2) = speed_function_table(j,1)
3545  speed_function_table(j,3) = speed_function_table(j,1)
3546 
3547  ENDDO
3548 
3549  CLOSE(11)
3550 
3551  ENDIF
3552 
3553  ENDIF
3554 
3555 ! for the ANISOTROPIC materials
3556  IF (material_behaviour .EQ. "ANISOTROPIC") THEN
3557 
3558 ! if it comes with FIXED format
3559  IF (input_type_for_speed_function .EQ. "FIXED") THEN
3560 
3561  DO i=1,total_number_of_nodes
3562  DO j=1,3
3563  speed_function_table(i,j) = speed_function_along_eigen_vector(j)
3564  ENDDO
3565  ENDDO
3566 
3567  ENDIF
3568 
3569 ! if it comes with FILE format
3570  IF (input_type_for_speed_function .EQ. "FILE") THEN
3571 
3572  text_length = index(input_file_name,' ') - 1
3573 
3574  string = input_file_name(1:text_length)//".cond"
3575 
3576  OPEN (11,file=string)
3577  READ (11,*) n
3578 
3579  DO i=1,n
3580 
3581  READ(11,*) j,speed_function_table(j,1),speed_function_table(j,2),&
3582  &speed_function_table(j,3)
3583 
3584  ENDDO
3585 
3586  CLOSE(11)
3587 
3588  ENDIF
3589 
3590  ENDIF
3591 
3592  ENDIF
3593 
3594  ! to set RAW_INDEX and COLUMN_INDEX
3595  raw_index(1) = 0
3596  DO i=1,total_number_of_nodes
3597 
3598  raw_index(i+1) = raw_index(i) + connectivity_number(i)
3599  DO j = 1,connectivity_number(i)
3600  column_index(raw_index(i)+j) = connectivity_list(i,j)
3601  ENDDO
3602 
3603  ENDDO
3604 
3605  DO i=1,total_number_of_nodes
3606 
3607  DO j = raw_index(i)+1,raw_index(i+1)
3608  DO k = 1,3
3609  speed_function_table_on_connectivity(j,k) = &
3610  & (speed_function_table(i,k)+speed_function_table(column_index(j),k))/2.0_dp
3611  ENDDO
3612  DO k = 1,9
3613  conductivity_tensor_on_connectivity(j,k) = &
3614  & (conductivity_tensor(i,k)+conductivity_tensor(column_index(j),k))/2.0_dp
3615  ENDDO
3616  ENDDO
3617 
3618  ENDDO
3619 
3620 ! EXITS("GENERATE_STATUS_MASK")
3621 ! RETURN
3622 999 errorsexits("GENERATE_STATUS_MASK",err,error)
3623  RETURN
3624 
3625  END SUBROUTINE pre_process_information
3626 
3627 
3628  !
3629  !================================================================================================================================
3630  !
3631 
3632 
3633  SUBROUTINE solve_problem_fmm(TOTAL_NUMBER_OF_NODES,NODE_LIST,CONDUCTIVITY_TENSOR,SPEED_FUNCTION_TABLE,&
3634  &seed_value,connectivity_number,connectivity_list,status_mask)
3635 
3636  !Argument variables
3637  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_NODES
3638 
3639  REAL(DP), ALLOCATABLE :: NODE_LIST(:,:)
3640  REAL(DP), ALLOCATABLE :: SPEED_FUNCTION_TABLE(:,:)
3641  REAL(DP), ALLOCATABLE :: CONDUCTIVITY_TENSOR(:,:)
3642  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_LIST(:,:)
3643 
3644  REAL(DP), ALLOCATABLE :: SEED_VALUE(:)
3645  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_NUMBER(:)
3646  CHARACTER (LEN=10), ALLOCATABLE :: STATUS_MASK(:)
3647 
3648  !Local Variables
3649  INTEGER(INTG) :: I,J
3650  INTEGER(INTG) :: LOOP_NUMBER, CHANGED_STATUS, MIN_TRIAL_NODE, TRIAL_STATUS
3651  REAL(DP), DIMENSION(3) :: DISTANCE_VECTOR,MV
3652  REAL(DP), DIMENSION(3,3) :: CONDUCTIVITY_MATRIX,F,FT,MFT,FMFT
3653  REAL(DP) :: MIN_TRIAL_VALUE,VMV,MINIMUM_DATA,TIME_ITER,TIME_NEW,CONDUCTION_RATIO
3654  INTEGER(INTG) :: Err
3655  TYPE(varying_string) :: Error
3656  err = 0
3657  !Start Program
3658 
3659  CALL generate_status_mask(total_number_of_nodes,seed_value,status_mask,err)
3660 
3661  min_trial_value = 1000
3662  DO i=1,total_number_of_nodes
3663 
3664  IF (status_mask(i) == "SEED POINT") THEN
3665 
3666  IF (seed_value(i) .lt. min_trial_value) THEN
3667  min_trial_value=seed_value(i)
3668  min_trial_node = i
3669  trial_status = 1
3670  ENDIF
3671 
3672  ENDIF
3673 
3674  ENDDO
3675 
3676 ! MAIN LOOP
3677  loop_number = 0
3678 
3679  DO WHILE (trial_status .eq. 1)
3680 
3681  trial_status = 0
3682  loop_number = loop_number + 1
3683  print *,"Running in loop number",loop_number
3684 
3685  ! CALL ASSIGN_MIN_TRIAL_TO_KNOWN(STATUS_MASK,MIN_TRIAL_NODE)
3686  status_mask(min_trial_node) = "KNOWN"
3687 
3688  DO i=1,connectivity_number(min_trial_node)
3689  time_new=1000
3690 
3691  DO j=1,connectivity_number(connectivity_list(min_trial_node,i))
3692  IF (status_mask(connectivity_list(connectivity_list(min_trial_node,i),j)) == "KNOWN") THEN
3693 
3694  distance_vector=(&
3695  &/node_list(connectivity_list(min_trial_node,i),1)&
3696  &-node_list(connectivity_list(connectivity_list(min_trial_node,i),j),1)&
3697  &,node_list(connectivity_list(min_trial_node,i),2)&
3698  &-node_list(connectivity_list(connectivity_list(min_trial_node,i),j),2)&
3699  &,node_list(connectivity_list(min_trial_node,i),3)&
3700  &-node_list(connectivity_list(connectivity_list(min_trial_node,i),j),3)/)
3701 
3702 ! CONDUCTION_RATIO=SPEED_FUNCTION_TABLE(CONNECTIVITY_LIST(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),J),1)/SPEED_FUNCTION_TABLE(CONNECTIVITY_LIST(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),J),2)
3703 ! CONDUCTIVITY_MATRIX=RESHAPE(SOURCE = (/1.0_DP,0.0_DP,0.0_DP,0.0_DP,CONDUCTION_RATIO*CONDUCTION_RATIO,0.0_DP,0.0_DP,0.0_DP,CONDUCTION_RATIO*CONDUCTION_RATIO/), SHAPE = (/3,3/))
3704  conduction_ratio= &
3705  &speed_function_table(connectivity_list(connectivity_list(min_trial_node,i),j),2)/&
3706  &speed_function_table(connectivity_list(connectivity_list(min_trial_node,i),j),1)
3707 ! &SPEED_FUNCTION_TABLE(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),CONNECTIVITY_LIST(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),J),2)
3708 
3709  conductivity_matrix=reshape(source=(/1.0_dp,0.0_dp,0.0_dp,&
3710  &0.0_dp,conduction_ratio*conduction_ratio,0.0_dp,&
3711  &0.0_dp,0.0_dp,conduction_ratio*conduction_ratio/), shape = (/3,3/))
3712 
3713 
3714 ! CALL LOAD_MATRIX((CONDUCTION_TENSOR(MIN_TRIAL_NODE,K),K=1,9),F)
3715  f=reshape(source =(/conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),1),&
3716  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),2),&
3717  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),3),&
3718  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),4),&
3719  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),5),&
3720  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),6),&
3721  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),7),&
3722  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),8),&
3723  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),9)/&
3724  &),shape = (/3,3/))
3725 
3726  CALL matrix_transpose(f,ft,err,error,*999)
3727  CALL matrix_product(conductivity_matrix,ft,mft,err,error,*999)
3728  CALL matrix_product(f,mft,fmft,err,error,*999)
3729 ! CALL INVERT(FMFT,INV_FMFT,DET,Err,Error,*999)
3730 
3731 ! PRINT *,F(1,1),F(1,2),F(1,3),F(2,1),F(2,2),F(2,3),F(3,1),F(3,2),F(3,3)
3732 
3733  CALL matrix_vector_product(fmft,distance_vector,mv,err,error,*999)
3734  CALL vector_vector_product(distance_vector,mv,vmv,err)
3735 
3736  time_iter=seed_value(connectivity_list(connectivity_list(min_trial_node,i),j))+sqrt(abs(vmv))*&
3737  &speed_function_table(connectivity_list(connectivity_list(min_trial_node,i),j),1)
3738 ! &SPEED_FUNCTION_TABLE(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),CONNECTIVITY_LIST(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),J),1)
3739 
3740 
3741  IF (time_iter .lt. time_new) THEN
3742  time_new = time_iter
3743  ENDIF
3744 
3745  ENDIF
3746  ENDDO
3747 
3748  IF (time_new .lt. seed_value(connectivity_list(min_trial_node,i))) THEN
3749  seed_value(connectivity_list(min_trial_node,i)) = time_new
3750  IF (status_mask(connectivity_list(min_trial_node,i)) .EQ. "KNOWN") THEN
3751  status_mask(connectivity_list(min_trial_node,i)) = "CHANGED"
3752  ENDIF
3753  IF (status_mask(connectivity_list(min_trial_node,i)) .EQ. "UNKNOWN") THEN
3754  status_mask(connectivity_list(min_trial_node,i)) = "SEED POINT"
3755  ENDIF
3756  ENDIF
3757 
3758  ENDDO
3759 
3760  minimum_data = 1000.0_dp
3761  changed_status = 0
3762 
3763  DO i=1,total_number_of_nodes
3764 
3765  IF (status_mask(i) .EQ. "CHANGED") THEN
3766  min_trial_node=i
3767  changed_status = 1
3768  trial_status = 1
3769  ENDIF
3770 
3771  IF (changed_status .EQ. 0.AND.status_mask(i) .EQ. "SEED POINT".AND.seed_value(i).LT.minimum_data) THEN
3772  minimum_data = seed_value(i)
3773  min_trial_node=i
3774  trial_status = 1
3775  ENDIF
3776 
3777  ENDDO
3778 
3779  ENDDO
3780 
3781 
3782 999 errorsexits("SOLVE_PROBLEM_FMM",err,error)
3783  RETURN
3784 
3785  END SUBROUTINE solve_problem_fmm
3786 
3787  !
3788  !================================================================================================================================
3789  !
3790 
3791 
3792  SUBROUTINE solve_problem_fmm_connectivity(TOTAL_NUMBER_OF_NODES,NODE_LIST,CONDUCTIVITY_TENSOR_ON_CONNECTIVITY,&
3793  &speed_function_table_on_connectivity,raw_index,column_index,total_number_of_connectivity,&
3794  &seed_value,status_mask)
3795 
3796  !Argument variables
3797  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_NODES
3798  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_CONNECTIVITY
3799 
3800  REAL(DP), ALLOCATABLE :: NODE_LIST(:,:)
3801  REAL(DP), ALLOCATABLE :: SPEED_FUNCTION_TABLE_ON_CONNECTIVITY(:,:)
3802  REAL(DP), ALLOCATABLE :: CONDUCTIVITY_TENSOR_ON_CONNECTIVITY(:,:)
3803  INTEGER(INTG), ALLOCATABLE :: COLUMN_INDEX(:)
3804  INTEGER(INTG), ALLOCATABLE :: RAW_INDEX(:)
3805  REAL(DP), ALLOCATABLE :: SEED_VALUE(:)
3806  CHARACTER (LEN=10), ALLOCATABLE :: STATUS_MASK(:)
3807 
3808  !Local Variables
3809  INTEGER(INTG) :: I,J
3810  INTEGER(INTG) :: LOOP_NUMBER, CHANGED_STATUS, MIN_TRIAL_NODE, TRIAL_STATUS
3811  REAL(DP), DIMENSION(3) :: DISTANCE_VECTOR,MV
3812  REAL(DP), DIMENSION(2) :: CONDUCTION_RATIO
3813  REAL(DP), DIMENSION(3,3) :: CONDUCTIVITY_MATRIX,F,FT,MFT,FMFT
3814  REAL(DP) :: MIN_TRIAL_VALUE,VMV,MINIMUM_DATA,TIME_ITER,TIME_NEW
3815  INTEGER(INTG) :: Err
3816  TYPE(varying_string) :: Error
3817  err = 0
3818  !Start Program
3819 
3820  CALL generate_status_mask(total_number_of_nodes,seed_value,status_mask,err)
3821 
3822  min_trial_value = 1000
3823  DO i=1,total_number_of_nodes
3824 
3825  IF (status_mask(i) == "SEED POINT") THEN
3826 
3827  IF (seed_value(i) .lt. min_trial_value) THEN
3828  min_trial_value=seed_value(i)
3829  min_trial_node = i
3830  trial_status = 1
3831  ENDIF
3832 
3833  ENDIF
3834 
3835  ENDDO
3836 
3837 ! ::::::: MAIN LOOP ::::::: P_NODE_NUMBER=CONNECTIVITY_LIST(MIN_TRIAL_NODE,I) PP_NODE_NUMBER=CONNECTIVITY_LIST(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),J)
3838  loop_number = 0
3839 
3840  DO WHILE (trial_status .eq. 1)
3841 
3842  trial_status = 0
3843  loop_number = loop_number + 1
3844  print *,"Running in loop number",loop_number
3845 
3846  ! CALL ASSIGN_MIN_TRIAL_TO_KNOWN(STATUS_MASK,MIN_TRIAL_NODE)
3847  status_mask(min_trial_node) = "KNOWN"
3848 
3849  DO i=raw_index(min_trial_node)+1,raw_index(min_trial_node+1)
3850  time_new=1000
3851 
3852  DO j=raw_index(column_index(i))+1,raw_index(column_index(i)+1)
3853  IF (status_mask(column_index(j)) == "KNOWN") THEN
3854 
3855  distance_vector=(/node_list(column_index(i),1)-node_list(column_index(j),1)&
3856  &,node_list(column_index(i),2)-node_list(column_index(j),2)&
3857  &,node_list(column_index(i),3)-node_list(column_index(j),3)/)
3858 
3859  conduction_ratio(1)= speed_function_table_on_connectivity(j,2)/speed_function_table_on_connectivity(j,1)
3860  conduction_ratio(2)= speed_function_table_on_connectivity(j,3)/speed_function_table_on_connectivity(j,1)
3861 
3862  conductivity_matrix=reshape(source=(/1.0_dp,0.0_dp,0.0_dp,&
3863  &0.0_dp,conduction_ratio(1)*conduction_ratio(1),0.0_dp,&
3864  &0.0_dp,0.0_dp,conduction_ratio(2)*conduction_ratio(2)/), shape = (/3,3/))
3865 
3866 
3867 ! CALL LOAD_MATRIX((CONDUCTION_TENSOR(MIN_TRIAL_NODE,K),K=1,9),F)
3868  f=reshape(source =(/conductivity_tensor_on_connectivity(j,1),&
3869  &conductivity_tensor_on_connectivity(j,2),&
3870  &conductivity_tensor_on_connectivity(j,3),&
3871  &conductivity_tensor_on_connectivity(j,4),&
3872  &conductivity_tensor_on_connectivity(j,5),&
3873  &conductivity_tensor_on_connectivity(j,6),&
3874  &conductivity_tensor_on_connectivity(j,7),&
3875  &conductivity_tensor_on_connectivity(j,8),&
3876  &conductivity_tensor_on_connectivity(j,9)/),shape = (/3,3/))
3877 
3878  CALL matrix_transpose(f,ft,err,error,*999)
3879  CALL matrix_product(conductivity_matrix,ft,mft,err,error,*999)
3880  CALL matrix_product(f,mft,fmft,err,error,*999)
3881 ! CALL INVERT(FMFT,INV_FMFT,DET,Err,Error,*999)
3882 
3883  CALL matrix_vector_product(fmft,distance_vector,mv,err,error,*999)
3884  CALL vector_vector_product(distance_vector,mv,vmv,err)
3885 
3886  time_iter=seed_value(column_index(j))+sqrt(abs(vmv))*speed_function_table_on_connectivity(j,1)
3887 
3888  IF (time_iter .lt. time_new) THEN
3889  time_new = time_iter
3890  ENDIF
3891 
3892  ENDIF
3893  ENDDO
3894 
3895  IF (time_new .LT. seed_value(column_index(i))) THEN
3896  seed_value(column_index(i)) = time_new
3897  IF (status_mask(column_index(i)) .EQ. "KNOWN") THEN
3898  status_mask(column_index(i)) = "CHANGED"
3899  ENDIF
3900  IF (status_mask(column_index(i)) .EQ. "UNKNOWN") THEN
3901  status_mask(column_index(i)) = "SEED POINT"
3902  ENDIF
3903  ENDIF
3904 
3905  ENDDO
3906 
3907  minimum_data = 1000.0_dp
3908  changed_status = 0
3909 
3910  DO i=1,total_number_of_nodes
3911 
3912  IF (status_mask(i) .EQ. "CHANGED") THEN
3913  min_trial_node=i
3914  changed_status = 1
3915  trial_status = 1
3916  ENDIF
3917 
3918  IF (changed_status .EQ. 0 .AND. status_mask(i) .EQ. "SEED POINT" .AND. seed_value(i) .LT. minimum_data) THEN
3919  minimum_data = seed_value(i)
3920  min_trial_node=i
3921  trial_status = 1
3922  ENDIF
3923 
3924  ENDDO
3925 
3926  ENDDO
3927 
3928 
3929 999 errorsexits("SOLVE_PROBLEM_FMM_CONNECTIVITY",err,error)
3930  RETURN
3931 
3932  END SUBROUTINE solve_problem_fmm_connectivity
3933 
3934 
3935  !
3936  !================================================================================================================================
3937  !
3938 
3939  SUBROUTINE solve_problem_geodesic_connectivity(TOTAL_NUMBER_OF_NODES,NODE_LIST,CONDUCTIVITY_TENSOR_ON_CONNECTIVITY,&
3940  &speed_function_table_on_connectivity,raw_index,column_index,total_number_of_connectivity,&
3941  &seed_value,status_mask,trace_node)
3942 
3943  !Argument variables
3944  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_NODES
3945  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_CONNECTIVITY
3946  REAL(DP), ALLOCATABLE :: NODE_LIST(:,:)
3947  REAL(DP), ALLOCATABLE :: SPEED_FUNCTION_TABLE_ON_CONNECTIVITY(:,:)
3948  REAL(DP), ALLOCATABLE :: CONDUCTIVITY_TENSOR_ON_CONNECTIVITY(:,:)
3949  INTEGER(INTG), ALLOCATABLE :: COLUMN_INDEX(:)
3950  INTEGER(INTG), ALLOCATABLE :: RAW_INDEX(:)
3951  REAL(DP), ALLOCATABLE :: SEED_VALUE(:)
3952  INTEGER(INTG), ALLOCATABLE :: TRACE_NODE(:)
3953  CHARACTER (LEN=10), ALLOCATABLE :: STATUS_MASK(:)
3954 
3955  !Local Variables
3956  REAL(DP), ALLOCATABLE :: CONNECTIVITY_WEIGHT(:)
3957  INTEGER(INTG) :: I,J
3958  INTEGER(INTG) :: CHANGED_STATUS,MIN_TRIAL_NODE,TRIAL_STATUS,TRACE_NODE_NUMBER
3959  REAL(DP), DIMENSION(3) :: DISTANCE_VECTOR,MV
3960  REAL(DP), DIMENSION(3,3) :: CONDUCTIVITY_MATRIX,F,FT,MFT,FMFT
3961  REAL(DP) :: MIN_TRIAL_VALUE,VMV,MINIMUM_DATA,TIME_ITER,TIME_NEW
3962  REAL(DP), DIMENSION(2) :: CONDUCTION_RATIO
3963  INTEGER(INTG) :: Err
3964  TYPE(varying_string) :: Error
3965  err = 0
3966 
3967  !initialize data
3968  DO i=1,total_number_of_nodes
3969  trace_node(i)=0
3970  ENDDO
3971 
3972  !Start Program
3973  ALLOCATE(connectivity_weight(total_number_of_connectivity),stat=err)
3974  DO i=1,total_number_of_nodes
3975  DO j=raw_index(i)+1,raw_index(i+1)
3976  distance_vector=(/node_list(i,1)-node_list(column_index(j),1)&
3977  &,node_list(i,2)-node_list(column_index(j),2)&
3978  &,node_list(i,3)-node_list(column_index(j),3)/)
3979 
3980  conduction_ratio(1) = speed_function_table_on_connectivity(j,2)/speed_function_table_on_connectivity(j,1)
3981  conduction_ratio(2) = speed_function_table_on_connectivity(j,3)/speed_function_table_on_connectivity(j,1)
3982 
3983  conductivity_matrix=reshape(source=(/1.0_dp,0.0_dp,0.0_dp,&
3984  &0.0_dp,conduction_ratio(1)*conduction_ratio(1),0.0_dp,&
3985  &0.0_dp,0.0_dp,conduction_ratio(2)*conduction_ratio(2)/), shape = (/3,3/))
3986 
3987  f=reshape(source =(/conductivity_tensor_on_connectivity(j,1),&
3988  &conductivity_tensor_on_connectivity(j,2),&
3989  &conductivity_tensor_on_connectivity(j,3),&
3990  &conductivity_tensor_on_connectivity(j,4),&
3991  &conductivity_tensor_on_connectivity(j,5),&
3992  &conductivity_tensor_on_connectivity(j,6),&
3993  &conductivity_tensor_on_connectivity(j,7),&
3994  &conductivity_tensor_on_connectivity(j,8),&
3995  &conductivity_tensor_on_connectivity(j,9)/),shape = (/3,3/))
3996 
3997  CALL matrix_transpose(f,ft,err,error,*999)
3998  CALL matrix_product(conductivity_matrix,ft,mft,err,error,*999)
3999  CALL matrix_product(f,mft,fmft,err,error,*999)
4000 
4001  CALL matrix_vector_product(fmft,distance_vector,mv,err,error,*999)
4002  CALL vector_vector_product(distance_vector,mv,vmv,err)
4003 
4004  connectivity_weight(j)=sqrt(abs(vmv))*speed_function_table_on_connectivity(j,1)
4005  ENDDO
4006  ENDDO
4007 
4008  CALL generate_status_mask(total_number_of_nodes,seed_value,status_mask,err)
4009 
4010  min_trial_value = 1000
4011  DO i=1,total_number_of_nodes
4012 
4013  IF (status_mask(i) == "SEED POINT") THEN
4014 
4015  IF (seed_value(i) .lt. min_trial_value) THEN
4016  min_trial_value=seed_value(i)
4017  min_trial_node = i
4018  trial_status = 1
4019  ENDIF
4020 
4021  ENDIF
4022 
4023  ENDDO
4024 
4025 ! ::::::: MAIN LOOP ::::::: P_NODE_NUMBER=CONNECTIVITY_LIST(MIN_TRIAL_NODE,I) PP_NODE_NUMBER=CONNECTIVITY_LIST(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),J)
4026 ! LOOP_NUMBER = 0
4027  print *,"Running GEODESIC Solver ...."
4028 
4029  DO WHILE (trial_status .eq. 1)
4030 
4031  trial_status = 0
4032  ! CALL ASSIGN_MIN_TRIAL_TO_KNOWN(STATUS_MASK,MIN_TRIAL_NODE)
4033  status_mask(min_trial_node) = "KNOWN"
4034 
4035  DO i=raw_index(min_trial_node)+1,raw_index(min_trial_node+1)
4036  time_new=1000
4037 
4038  DO j=raw_index(column_index(i))+1,raw_index(column_index(i)+1)
4039  IF (status_mask(column_index(j)) == "KNOWN") THEN
4040 
4041  time_iter=seed_value(column_index(j))+connectivity_weight(j)
4042 
4043  IF (time_iter .lt. time_new) THEN
4044  time_new = time_iter
4045  trace_node_number=column_index(j)
4046  ENDIF
4047 
4048  ENDIF
4049  ENDDO
4050 
4051  IF (time_new .LT. seed_value(column_index(i))) THEN
4052  seed_value(column_index(i)) = time_new
4053  trace_node(column_index(i)) = trace_node_number
4054  IF (status_mask(column_index(i)) .EQ. "KNOWN") THEN
4055  status_mask(column_index(i)) = "CHANGED"
4056  ENDIF
4057  IF (status_mask(column_index(i)) .EQ. "UNKNOWN") THEN
4058  status_mask(column_index(i)) = "SEED POINT"
4059  ENDIF
4060  ENDIF
4061 
4062  ENDDO
4063 
4064  minimum_data = 1000.0_dp
4065  changed_status = 0
4066 
4067  DO i=1,total_number_of_nodes
4068 
4069  IF (status_mask(i) .EQ. "CHANGED") THEN
4070  min_trial_node = i
4071  changed_status = 1
4072  trial_status = 1
4073  ENDIF
4074 
4075  IF (changed_status .EQ. 0 .AND. status_mask(i) .EQ. "SEED POINT" .AND. seed_value(i) .LT. minimum_data) THEN
4076  minimum_data = seed_value(i)
4077  min_trial_node = i
4078  trial_status = 1
4079  ENDIF
4080 
4081  ENDDO
4082 
4083  ENDDO
4084 
4085 999 errorsexits("SOLVE_PROBLEM_GEODESIC_CONNECTIVITY",err,error)
4086  RETURN
4087 
4088  END SUBROUTINE solve_problem_geodesic_connectivity
4089 
4090  !
4091  !================================================================================================================================
4092  !
4093 
4094 
4095  SUBROUTINE solve_problem_geodesic(TOTAL_NUMBER_OF_NODES,NODE_LIST,CONDUCTIVITY_TENSOR,SPEED_FUNCTION_TABLE,&
4096  & seed_value,connectivity_number,connectivity_list,status_mask,trace_node)
4097 
4098  !Argument variables
4099  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_NODES
4100 
4101  REAL(DP), ALLOCATABLE :: NODE_LIST(:,:)
4102  REAL(DP), ALLOCATABLE :: SPEED_FUNCTION_TABLE(:,:)
4103  REAL(DP), ALLOCATABLE :: CONDUCTIVITY_TENSOR(:,:)
4104  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_LIST(:,:)
4105 
4106  REAL(DP), ALLOCATABLE :: SEED_VALUE(:)
4107  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_NUMBER(:)
4108  INTEGER(INTG), ALLOCATABLE :: TRACE_NODE(:)
4109  CHARACTER (LEN=10), ALLOCATABLE :: STATUS_MASK(:)
4110 
4111  !Local Variables
4112  INTEGER(INTG) :: I,J
4113  INTEGER(INTG) :: LOOP_NUMBER,CHANGED_STATUS,MIN_TRIAL_NODE,TRIAL_STATUS,TRACE_NODE_NUMBER
4114  REAL(DP), DIMENSION(3) :: DISTANCE_VECTOR,MV
4115  REAL(DP), DIMENSION(3,3) :: CONDUCTIVITY_MATRIX,F,FT,MFT,FMFT
4116  REAL(DP) :: MIN_TRIAL_VALUE,VMV,MINIMUM_DATA,TIME_ITER,TIME_NEW,CONDUCTION_RATIO
4117  INTEGER(INTG) :: Err
4118  TYPE(varying_string) :: Error
4119  err = 0
4120 
4121  !Start Program
4122 
4123  CALL generate_status_mask(total_number_of_nodes,seed_value,status_mask,err)
4124 
4125  min_trial_value = 1000
4126  DO i=1,total_number_of_nodes
4127 
4128  IF (status_mask(i) == "SEED POINT") THEN
4129 
4130  IF (seed_value(i) .lt. min_trial_value) THEN
4131  min_trial_value=seed_value(i)
4132  min_trial_node = i
4133  trial_status = 1
4134  ENDIF
4135 
4136  ENDIF
4137 
4138  ENDDO
4139 
4140 ! 4444444444444444 MAIN LOOP 55555555555555555 P_NODE_NUMBER=CONNECTIVITY_LIST(MIN_TRIAL_NODE,I) PP_NODE_NUMBER=CONNECTIVITY_LIST(CONNECTIVITY_LIST(MIN_TRIAL_NODE,I),J)
4141  loop_number = 0
4142 
4143  DO WHILE (trial_status .eq. 1)
4144 
4145  trial_status = 0
4146  loop_number = loop_number + 1
4147  print *,"Running in loop number",loop_number
4148 
4149  ! CALL ASSIGN_MIN_TRIAL_TO_KNOWN(STATUS_MASK,MIN_TRIAL_NODE)
4150  status_mask(min_trial_node) = "KNOWN"
4151 
4152  DO i=1,connectivity_number(min_trial_node)
4153  time_new=1000
4154 
4155  DO j=1,connectivity_number(connectivity_list(min_trial_node,i))
4156  IF (status_mask(connectivity_list(connectivity_list(min_trial_node,i),j)) == "KNOWN") THEN
4157 
4158  distance_vector=(&
4159  &/node_list(connectivity_list(min_trial_node,i),1)&
4160  &-node_list(connectivity_list(connectivity_list(min_trial_node,i),j),1)&
4161  &,node_list(connectivity_list(min_trial_node,i),2)&
4162  &-node_list(connectivity_list(connectivity_list(min_trial_node,i),j),2)&
4163  &,node_list(connectivity_list(min_trial_node,i),3)&
4164  &-node_list(connectivity_list(connectivity_list(min_trial_node,i),j),3)/)
4165 
4166  conduction_ratio=speed_function_table(connectivity_list(connectivity_list(min_trial_node,i),j),2)/&
4167  &speed_function_table(connectivity_list(connectivity_list(min_trial_node,i),j),1)
4168 
4169  conductivity_matrix=reshape(source=(/1.0_dp,0.0_dp,0.0_dp,&
4170  &0.0_dp,conduction_ratio*conduction_ratio,0.0_dp,&
4171  &0.0_dp,0.0_dp,conduction_ratio*conduction_ratio/), shape = (/3,3/))
4172 
4173 ! CALL LOAD_MATRIX((CONDUCTION_TENSOR(MIN_TRIAL_NODE,K),K=1,9),F)
4174  f=reshape(source =(/conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),1),&
4175  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),2),&
4176  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),3),&
4177  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),4),&
4178  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),5),&
4179  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),6),&
4180  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),7),&
4181  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),8),&
4182  &conductivity_tensor(connectivity_list(connectivity_list(min_trial_node,i),j),9)/&
4183  &),shape = (/3,3/))
4184 
4185  CALL matrix_transpose(f,ft,err,error,*999)
4186  CALL matrix_product(conductivity_matrix,ft,mft,err,error,*999)
4187  CALL matrix_product(f,mft,fmft,err,error,*999)
4188 
4189  CALL matrix_vector_product(fmft,distance_vector,mv,err,error,*999)
4190  CALL vector_vector_product(distance_vector,mv,vmv,err)
4191 
4192  time_iter=seed_value(connectivity_list(connectivity_list(min_trial_node,i),j))+&
4193  &sqrt(abs(vmv))*speed_function_table(connectivity_list(connectivity_list(min_trial_node,i),j),1)
4194 
4195  IF (time_iter .lt. time_new) THEN
4196  time_new = time_iter
4197  trace_node_number=connectivity_list(connectivity_list(min_trial_node,i),j)
4198  ENDIF
4199 
4200  ENDIF
4201  ENDDO
4202 
4203  IF (time_new .lt. seed_value(connectivity_list(min_trial_node,i))) THEN
4204  seed_value(connectivity_list(min_trial_node,i)) = time_new
4205  trace_node(connectivity_list(min_trial_node,i))=trace_node_number
4206  IF (status_mask(connectivity_list(min_trial_node,i)) .EQ. "KNOWN") THEN
4207  status_mask(connectivity_list(min_trial_node,i)) = "CHANGED"
4208  ENDIF
4209  IF (status_mask(connectivity_list(min_trial_node,i)) .EQ. "UNKNOWN") THEN
4210  status_mask(connectivity_list(min_trial_node,i)) = "SEED POINT"
4211  ENDIF
4212  ENDIF
4213 
4214  ENDDO
4215 
4216  minimum_data = 1000
4217  changed_status = 0
4218 
4219  DO i=1,total_number_of_nodes
4220 
4221  IF (status_mask(i) .EQ. "CHANGED") THEN
4222  min_trial_node=i
4223  changed_status = 1
4224  trial_status = 1
4225  ENDIF
4226 
4227  IF (changed_status .EQ. 0.AND.status_mask(i) .EQ. "SEED POINT".AND.seed_value(i).LT.minimum_data) THEN
4228  minimum_data = seed_value(i)
4229  min_trial_node=i
4230  trial_status = 1
4231  ENDIF
4232 
4233  ENDDO
4234 
4235  ENDDO
4236 
4237 
4238 999 errorsexits("SOLVE_PROBLEM_GEODESIC",err,error)
4239  RETURN
4240 
4241  END SUBROUTINE solve_problem_geodesic
4242 
4243 
4244  !
4245  !================================================================================================================================
4246  !
4247 
4249  SUBROUTINE generate_status_mask(TOTAL_NUMBER_OF_NODES,SEED_VALUE,STATUS_MASK,Err)
4250 
4251  !subroutine variables
4252  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_NODES
4253  REAL(DP), ALLOCATABLE :: SEED_VALUE(:)
4254  CHARACTER (LEN=10), ALLOCATABLE :: STATUS_MASK(:)
4255 
4256  !Argument variables
4257  INTEGER(INTG) :: Err
4258  ! TYPE(VARYING_STRING) :: LOCAL_ERROR !<The error string
4259 
4260  !Local variables
4261  INTEGER(INTG) :: I
4262 
4263 ! ENTERS("GENERATE_STATUS_MASK",Err,Error,*999)
4264 
4265  DO i=1,total_number_of_nodes
4266 
4267  IF (seed_value(i) .LT. 100.0_dp) THEN
4268  status_mask(i) = "SEED POINT"
4269  ELSE
4270  status_mask(i) = "UNKNOWN"
4271  ENDIF
4272 
4273  ENDDO
4274 
4275 ! EXITS("GENERATE_STATUS_MASK")
4276 ! RETURN
4277 !999 ERRORSEXITS("GENERATE_STATUS_MASK",ERR,ERROR)
4278 ! RETURN 1
4279 
4280  END SUBROUTINE generate_status_mask
4281 
4282 
4283  !
4284  !================================================================================================================================
4285  !
4286 
4288  SUBROUTINE find_minimax(A,N,MIN_VALUE,MAX_VALUE,Err)
4289 
4290  !Argument variables
4291  REAL(DP), ALLOCATABLE :: A(:)
4292 
4293  REAL(DP), INTENT(OUT) :: MIN_VALUE
4294  REAL(DP), INTENT(OUT) :: MAX_VALUE
4295  INTEGER(INTG), INTENT(IN) :: N
4296  INTEGER(INTG) :: ERR
4297  ! TYPE(VARYING_STRING) :: LOCAL_ERROR !<The error string
4298  !Local variables
4299  INTEGER(INTG) :: I
4300 
4301 ! ENTERS("FIND_MINIMAX",Err,Error,*999)
4302 
4303  IF(SIZE(a,1).GT.2) THEN
4304  min_value = a(1)
4305  max_value = a(1)
4306  DO i=2,n
4307  IF (min_value .GT. a(i)) THEN
4308  min_value = a(i)
4309  ENDIF
4310  IF (max_value .LT. a(i)) THEN
4311  max_value = a(i)
4312  ENDIF
4313  ENDDO
4314 ! ELSE
4315 ! CALL FlagError("Invalid matrix sizes.",Err)
4316  ENDIF
4317 
4318 ! EXITS("FIND_MINIMAX")
4319 ! RETURN
4320 !999 ERRORSEXITS("FIND_MINIMAX",ERR,ERROR)
4321 ! RETURN 1
4322 
4323  END SUBROUTINE find_minimax
4324 
4325  !
4326  !================================================================================================================================
4327  !
4328 
4329 
4331  SUBROUTINE vector_vector_product(A,B,C,Err)
4332 
4333  !Argument variables
4334  REAL(DP), INTENT(IN) :: A(3)
4335  REAL(DP), INTENT(IN) :: B(3)
4336  REAL(DP), INTENT(OUT) :: C
4337  INTEGER(INTG) :: ERR
4338  ! TYPE(VARYING_STRING) :: LOCAL_ERROR !<The error string
4339  !Local variables
4340 
4341 ! ENTERS("VECTOR_VECTOR_PRODUCT",Err,Error,*999)
4342 
4343  IF(SIZE(a,1)==SIZE(b,1)) THEN
4344  SELECT CASE(SIZE(a,1))
4345  CASE(1)
4346  c=a(1)*b(1)
4347  CASE(2)
4348  c=a(1)*b(1)+a(2)*b(2)
4349  CASE(3)
4350  c=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
4351 ! CASE DEFAULT
4352 ! CALL FlagError("Invalid matrix size.",Err)
4353  END SELECT
4354 ! ELSE
4355 ! CALL FlagError("Invalid matrix sizes.",Err)
4356  ENDIF
4357 
4358 ! EXITS("VECTOR_VECTOR_PRODUCT")
4359 ! RETURN
4360 !999 ERRORSEXITS("VECTOR_VECTOR_PRODUCT",ERR,ERROR)
4361 ! RETURN 1
4362 
4363  END SUBROUTINE vector_vector_product
4364 
4365  !
4366  !================================================================================================================================
4367  !
4368 
4369 
4371  SUBROUTINE post_process_data(MATERIAL_BEHAVIOUR,OUTPUT_FILE_NAME,OUTPUT_FILE_FORMAT,TOTAL_NUMBER_OF_NODES,NODE_LIST,&
4372 &conductivity_tensor,speed_function_table,seed_value,connectivity_number,output_file_field_title,&
4373 &connectivity_list,element_list,total_number_of_elements,number_of_nodes_per_element,err)
4374 
4375  !subroutine variables
4376  REAL(DP), ALLOCATABLE :: NODE_LIST(:,:)
4377  REAL(DP), ALLOCATABLE :: SPEED_FUNCTION_TABLE(:,:)
4378  REAL(DP), ALLOCATABLE :: CONDUCTIVITY_TENSOR(:,:)
4379  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_LIST(:,:)
4380  INTEGER(INTG), ALLOCATABLE :: ELEMENT_LIST(:,:)
4381  REAL(DP), ALLOCATABLE :: SEED_VALUE(:)
4382  INTEGER(INTG), ALLOCATABLE :: CONNECTIVITY_NUMBER(:)
4383 
4384  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_ELEMENTS
4385  INTEGER(INTG), INTENT(IN) :: NUMBER_OF_NODES_PER_ELEMENT
4386  INTEGER(INTG), INTENT(IN) :: TOTAL_NUMBER_OF_NODES
4387  CHARACTER (LEN=300) :: OUTPUT_FILE_NAME
4388  CHARACTER (LEN=300) :: OUTPUT_FILE_FIELD_TITLE
4389  CHARACTER (LEN=10) :: OUTPUT_FILE_FORMAT
4390  CHARACTER (LEN=12) :: MATERIAL_BEHAVIOUR
4391  INTEGER(INTG) :: Err
4392 
4393  !Local variables
4394  CHARACTER (LEN=300) :: STRING
4395  INTEGER(INTG) :: TEXT_LENGTH
4396  INTEGER(INTG) :: I, J
4397 
4398 ! ENTERS("GENERATE_STATUS_MASK",Err,Error,*999)
4399 
4400 ! IF (OUTPUT_FILE_FORMAT .NE. "TABC") THEN
4401 ! CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,">> OUTPUT FILE FORMAT ERROR: identify a correct output format",Err)
4402 ! ENDIF
4403 
4404 
4405 ! in the case the OUTPUT is in TABC format (2)
4406  IF (output_file_format .EQ. "TABC") THEN
4407 
4408 ! EXPORT NODE TABLE list
4409  text_length = index(output_file_name,' ') - 1
4410  string = output_file_name(1:text_length)//".tabc"
4411 ! INQUIRE(FILE=STRING, EXIST=ex)
4412 ! PRINT *, STRING
4413  OPEN (12,file=string)
4414 ! OPEN (12,FILE=OUTPUT_FILE_NAME)
4415 
4416  WRITE(12,*)"VARIABLES=""NODE"",""X"",""Y"",""Z"",""U"",""V"",""W"",""Speed function along fibers"", &
4417  & ""Speed function in transverse direction"",""Time"""
4418  WRITE(12,*)"zone i=",total_number_of_nodes," , DATAPACKING=POINT"
4419 
4420  DO i=1,total_number_of_nodes
4421  WRITE(12,*) i,node_list(i,1),node_list(i,2),node_list(i,3),conductivity_tensor(i,1),conductivity_tensor(i,2),&
4422  &conductivity_tensor(i,3),speed_function_table(i,1),&
4423  &speed_function_table(i,2),speed_function_table(i,3),&
4424  &seed_value(i)
4425  ENDDO
4426 ! EXPORT NODE CONNECTIVITY list
4427  WRITE(12,*) "Connectivity"
4428  DO i=1,total_number_of_nodes
4429  WRITE(12,*) i,connectivity_number(i),(connectivity_list(i,j),j=1,connectivity_number(i))
4430 ! PRINT *,STRING,CONNECTIVITY_NUMBER(I),(CONNECTIVITY_LIST(I,J),J=1,CONNECTIVITY_NUMBER(I))
4431  ENDDO
4432  CLOSE(12)
4433 
4434  ENDIF
4435 
4436 ! in the case the OUTPUT is in VTK TETrahedral format (2)
4437  IF (output_file_format .EQ. "VTKTET") THEN
4438 
4439 ! rename to VTK and OPEN the file
4440  text_length = index(output_file_name,' ') - 1
4441  string = output_file_name(1:text_length)//".vtk"
4442  OPEN (12,file=string)
4443 
4444 ! export HEADER list
4445  WRITE(12,'(A)')"# vtk DataFile Version 3.0"
4446  WRITE(12,'(A)')"vtk output"
4447  WRITE(12,'(A)')"ASCII"
4448  WRITE(12,'(A)')"DATASET UNSTRUCTURED_GRID"
4449  WRITE(12,'(A,I8,A6)')"POINTS",total_number_of_nodes,"float"
4450 
4451 ! export NODAL POSITION list
4452  DO i=1,total_number_of_nodes
4453  WRITE(12,*) (node_list(i,j),j=1,3)
4454  ENDDO
4455 
4456 ! export ELEMENT list
4457  WRITE(12,'(A,I8,I8)')"CELLS",total_number_of_elements,total_number_of_elements*(number_of_nodes_per_element+1)
4458  DO i=1,total_number_of_elements
4459  WRITE(12,*) number_of_nodes_per_element,((element_list(i,j)-1),j=1,number_of_nodes_per_element)
4460  ENDDO
4461 
4462 ! export ELEMENT TYPE list
4463  WRITE(12,'(A,I8)')"CELL_TYPES",total_number_of_elements
4464  DO i=1,total_number_of_elements
4465  WRITE(12,'(A)') "10"
4466  ENDDO
4467 
4468 ! export CELL and POINT DATA list
4469  WRITE(12,'(A,I8)')"CELL_DATA",total_number_of_elements
4470  WRITE(12,'(A,I8)')"POINT_DATA",total_number_of_nodes
4471 
4472 ! export FIELD information
4473  WRITE(12,'(A,A)')"FIELD number"," 1"
4474  WRITE(12,'(A,I3,I8,A6)')output_file_field_title,1,total_number_of_nodes,"float"
4475  DO i=1,total_number_of_nodes
4476  WRITE(12,'(F15.10)') seed_value(i)
4477  ENDDO
4478 
4479 ! export VECTORS information
4480  WRITE(12,'(A,A,A6)') "VECTORS ","fiber_vector","float"
4481  DO i=1,total_number_of_nodes
4482  WRITE(12,'(3F8.5)') (conductivity_tensor(i,j),j=1,3)
4483  ENDDO
4484 
4485  CLOSE(12)
4486 
4487  ENDIF
4488 
4489 ! FILE="cmgui"
4490 ! METHOD="FORTRAN"
4491 
4492 ! EXPORT_FIELD=.TRUE.
4493 ! IF(EXPORT_FIELD) THEN
4494 ! WRITE(*,*)'Now export fields...'
4495 ! CALL FLUID_MECHANICS_IO_WRITE_CMGUI(REGION,FILE,Err)
4496 ! CALL FIELD_IO_NODES_EXPORT(REGION%FIELDS, FILE, METHOD, Err)
4497 ! CALL FIELD_IO_ELEMENTS_EXPORT(REGION%FIELDS, FILE, METHOD, Err)
4498 ! WRITE(*,*)'All fields exported...'
4499 ! ENDIF
4500 
4501 ! EXITS("GENERATE_STATUS_MASK")
4502 ! RETURN
4503 !999 ERRORSEXITS("GENERATE_STATUS_MASK",ERR,ERROR)
4504 ! RETURN 1
4505 
4506  END SUBROUTINE post_process_data
4507 
4508 
4509 END MODULE hamilton_jacobi_equations_routines
equations_set_constants::equations_set_setup_dependent_type
integer(intg), parameter equations_set_setup_dependent_type
Dependent variables.
Definition: equations_set_constants.f90:324
equations_set_constants::equations_set_fem_solution_method
integer(intg), parameter equations_set_fem_solution_method
Finite Element Method solution method.
Definition: equations_set_constants.f90:348
basis_routines
This module contains all basis function routines.
Definition: basis_routines.f90:45
equations_set_constants::equations_set_setup_materials_type
integer(intg), parameter equations_set_setup_materials_type
Materials setup.
Definition: equations_set_constants.f90:326
types::boundary_conditions_type
Contains information on the boundary conditions for the solver equations.
Definition: types.f90:1780
base_routines::enters
subroutine, public enters(NAME, ERR, ERROR,)
Records the entry into the named procedure and initialises the error code.
Definition: base_routines.f90:363
solver_routines::solvers_create_finish
subroutine, public solvers_create_finish(SOLVERS, ERR, ERROR,)
Finish the creation of solvers.
Definition: solver_routines.f90:22598
equations_set_constants::equations_set_hj_equation_three_dim_1
integer(intg), parameter equations_set_hj_equation_three_dim_1
u=x**2-2*y**2+z**2
Definition: equations_set_constants.f90:431
types::equations_mapping_type
Contains information on the equations mapping i.e., how field variable DOFS are mapped to the rows an...
Definition: types.f90:1681
types::equations_type
Contains information about the equations in an equations set.
Definition: types.f90:1735
equations_set_constants::equations_set_gfem_solution_method
integer(intg), parameter equations_set_gfem_solution_method
Grid-based Finite Element Method solution method.
Definition: equations_set_constants.f90:352
constants::global_deriv_s1_s2_s3
integer(intg), parameter global_deriv_s1_s2_s3
Cross derivative in the s1, s2 and s3 direction i.e., d^3u/ds1ds2ds3.
Definition: constants.f90:220
equations_set_constants::equations_set_generalised_hj_subtype
integer(intg), parameter equations_set_generalised_hj_subtype
Definition: equations_set_constants.f90:218
problem_constants::problem_setup_control_type
integer(intg), parameter problem_setup_control_type
Solver setup for a problem.
Definition: problem_constants.f90:263
problem_constants
This module handles all problem wide constants.
Definition: problem_constants.f90:45
control_loop_routines::control_loop_node
integer(intg), parameter, public control_loop_node
The identifier for a each "leaf" node in a control loop.
Definition: control_loop_routines.f90:67
maths::matrix_transpose
Returns the transpose of a matrix A in A^T.
Definition: maths.f90:191
constants::no_global_deriv
integer(intg), parameter no_global_deriv
No global derivative i.e., u.
Definition: constants.f90:213
strings::number_to_vstring
Converts a number to its equivalent varying string representation.
Definition: strings.f90:161
equations_routines::equations_create_start
subroutine, public equations_create_start(EQUATIONS_SET, EQUATIONS, ERR, ERROR,)
Start the creation of equations for the equation set.
Definition: equations_routines.f90:192
types::decomposition_type
Contains information on the mesh decomposition.
Definition: types.f90:1063
equations_matrices_routines::equations_matrices_create_start
subroutine, public equations_matrices_create_start(EQUATIONS, EQUATIONS_MATRICES, ERR, ERROR,)
Starts the creation of the equations matrices and rhs for the the equations.
Definition: equations_matrices_routines.f90:615
types::solver_type
Contains information on the type of solver to be used.
Definition: types.f90:2777
solver_routines::solver_petsc_library
integer(intg), parameter, public solver_petsc_library
PETSc solver library.
Definition: solver_routines.f90:104
solver_routines::solvers_number_set
subroutine, public solvers_number_set(SOLVERS, NUMBER_OF_SOLVERS, ERR, ERROR,)
Sets/changes the number of solvers.
Definition: solver_routines.f90:22805
hamilton_jacobi_equations_routines
This module handles all Hamilton-Jacobi equations routines.
Definition: Hamilton_Jacobi_equations_routines.f90:45
equations_matrices_routines
This module handles all equations matrix and rhs routines.
Definition: equations_matrices_routines.f90:45
solver_routines::solver_type_set
subroutine, public solver_type_set(SOLVER, SOLVE_TYPE, ERR, ERROR,)
Sets/changes the type for a solver.
Definition: solver_routines.f90:21345
equations_set_constants::equations_static
integer(intg), parameter equations_static
The equations are static and have no time dependence.
Definition: equations_set_constants.f90:376
types::equations_set_type
Contains information on an equations set.
Definition: types.f90:1941
equations_routines
This module handles all equations routines.
Definition: equations_routines.f90:45
equations_set_constants::equations_set_setup_source_type
integer(intg), parameter equations_set_setup_source_type
Source setup.
Definition: equations_set_constants.f90:327
strings
This module contains all string manipulation and transformation routines.
Definition: strings.f90:45
solver_routines::solvers_create_start
subroutine, public solvers_create_start(CONTROL_LOOP, SOLVERS, ERR, ERROR,)
Start the creation of a solvers for the control loop.
Definition: solver_routines.f90:22651
types::solvers_type
Contains information on the solvers to be used in a control loop.
Definition: types.f90:2805
constants::first_part_deriv
integer(intg), parameter first_part_deriv
First partial derivative i.e., du/ds.
Definition: constants.f90:178
equations_set_constants::equations_set_hj_equation_three_dim_2
integer(intg), parameter equations_set_hj_equation_three_dim_2
u=cos(x)*cosh(y)*z
Definition: equations_set_constants.f90:432
timer
This module contains routines for timing the program.
Definition: timer_f.f90:45
problem_constants::solver_equations_static
integer(intg), parameter solver_equations_static
Solver equations are static.
Definition: problem_constants.f90:309
equations_routines::equations_time_dependence_type_set
subroutine, public equations_time_dependence_type_set(EQUATIONS, TIME_DEPENDENCE_TYPE, ERR, ERROR,)
Sets/changes the time dependence type for equations.
Definition: equations_routines.f90:851
solver_routines::solver_equations_sparsity_type_set
subroutine, public solver_equations_sparsity_type_set(SOLVER_EQUATIONS, SPARSITY_TYPE, ERR, ERROR,)
Sets/changes the sparsity type for solver equations.
Definition: solver_routines.f90:7734
maths
This module contains all mathematics support routines.
Definition: maths.f90:45
solver_routines::solvers_solver_get
subroutine, public solvers_solver_get(SOLVERS, SOLVER_INDEX, SOLVER, ERR, ERROR,)
Returns a pointer to the specified solver in the list of solvers.
Definition: solver_routines.f90:22876
problem_constants::problem_hj_equation_type
integer(intg), parameter problem_hj_equation_type
Definition: problem_constants.f90:94
types::field_type
Contains information for a field defined on a region.
Definition: types.f90:1346
equations_matrices_routines::equations_matrices_full_matrices
integer(intg), parameter, public equations_matrices_full_matrices
Use fully populated equation matrices.
Definition: equations_matrices_routines.f90:91
equations_mapping_routines::equations_mapping_rhs_variable_type_set
subroutine, public equations_mapping_rhs_variable_type_set(EQUATIONS_MAPPING, RHS_VARIABLE_TYPE, ERR, ERROR,)
Sets the mapping between a dependent field variable and the equations set rhs vector.
Definition: equations_mapping_routines.f90:3170
problem_constants::solver_equations_linear
integer(intg), parameter solver_equations_linear
Solver equations are linear.
Definition: problem_constants.f90:301
constants::global_deriv_s2
integer(intg), parameter global_deriv_s2
First global derivative in the s2 direction i.e., du/ds2.
Definition: constants.f90:215
types::control_loop_type
Contains information on a control loop.
Definition: types.f90:3185
iso_varying_string
This module provides an iso_varying_string module, conformant to the API specified in ISO/IEC 1539-2:...
Definition: iso_varying_string.f90:62
solver_routines::solver_equations_create_finish
subroutine, public solver_equations_create_finish(SOLVER_EQUATIONS, ERR, ERROR,)
Finishes the process of creating solver equations.
Definition: solver_routines.f90:7166
solver_routines::solver_sparse_matrices
integer(intg), parameter, public solver_sparse_matrices
Use sparse solver matrices.
Definition: solver_routines.f90:360
constants::global_deriv_s2_s3
integer(intg), parameter global_deriv_s2_s3
Global Cross derivative in the s2 and s3 direction i.e., d^2u/ds2ds3.
Definition: constants.f90:219
solver_routines::solver_equations_create_start
subroutine, public solver_equations_create_start(SOLVER, SOLVER_EQUATIONS, ERR, ERROR,)
Starts the process of creating solver equations.
Definition: solver_routines.f90:7209
equations_set_constants::equations_set_hj_equation_two_dim_1
integer(intg), parameter equations_set_hj_equation_two_dim_1
u=x**2+y**2
Definition: equations_set_constants.f90:429
basis_routines::basis_default_quadrature_scheme
integer(intg), parameter, public basis_default_quadrature_scheme
Identifier for the default quadrature scheme.
Definition: basis_routines.f90:125
problem_constants::problem_setup_solvers_type
integer(intg), parameter problem_setup_solvers_type
Solver setup for a problem.
Definition: problem_constants.f90:264
equations_set_constants::equations_set_setup_equations_type
integer(intg), parameter equations_set_setup_equations_type
Equations setup.
Definition: equations_set_constants.f90:330
equations_set_constants::equations_set_standard_hj_subtype
integer(intg), parameter equations_set_standard_hj_subtype
Definition: equations_set_constants.f90:217
iso_varying_string::varying_string
Definition: iso_varying_string.f90:75
constants
This module contains all program wide constants.
Definition: constants.f90:45
solver_routines::solver_library_type_set
subroutine, public solver_library_type_set(SOLVER, SOLVER_LIBRARY_TYPE, ERR, ERROR,)
Sets/changes the type of library type to use for the solver.
Definition: solver_routines.f90:9067
maths::cross_product
Calculates the vector cross product of two vectors.
Definition: maths.f90:66
equations_mapping_routines::equationsmapping_linearmatricesnumberset
subroutine, public equationsmapping_linearmatricesnumberset(EQUATIONS_MAPPING, NUMBER_OF_LINEAR_EQUATIONS_MATRICES, ERR, ERROR,)
Sets/changes the number of linear equations matrices.
Definition: equations_mapping_routines.f90:2484
problem_constants::problem_setup_initial_type
integer(intg), parameter problem_setup_initial_type
Initial setup for a problem.
Definition: problem_constants.f90:262
equations_set_constants::equations_set_constant_source_poisson_subtype
integer(intg), parameter equations_set_constant_source_poisson_subtype
Definition: equations_set_constants.f90:221
equations_mapping_routines::equationsmapping_linearmatricesvariabletypesset
subroutine, public equationsmapping_linearmatricesvariabletypesset(EQUATIONS_MAPPING, LINEAR_MATRIX_VARIABLE_TYPES, ERR, ERROR,)
Sets the mapping between the dependent field variable types and the linear equations matrices...
Definition: equations_mapping_routines.f90:2716
types::problem_setup_type
Definition: types.f90:3215
solver_routines::solver_equations_linearity_type_set
subroutine, public solver_equations_linearity_type_set(SOLVER_EQUATIONS, LINEARITY_TYPE, ERR, ERROR,)
Sets/changes the linearity type for solver equations.
Definition: solver_routines.f90:7577
equations_set_constants::equations_set_setup_start_action
integer(intg), parameter equations_set_setup_start_action
Start setup action.
Definition: equations_set_constants.f90:338
problem_constants::problem_classical_field_class
integer(intg), parameter problem_classical_field_class
Definition: problem_constants.f90:58
base_routines::exits
subroutine, public exits(NAME)
Records the exit out of the named procedure.
Definition: base_routines.f90:504
control_loop_routines::control_loop_solvers_get
recursive subroutine, public control_loop_solvers_get(CONTROL_LOOP, SOLVERS, ERR, ERROR,)
Returns a pointer to the solvers for a control loop.
Definition: control_loop_routines.f90:1333
types
This module contains all type definitions in order to avoid cyclic module references.
Definition: types.f90:70
types::equations_matrices_type
Contains information on the equations matrices and vectors.
Definition: types.f90:1520
equations_matrices_routines::equations_matrix_fem_structure
integer(intg), parameter, public equations_matrix_fem_structure
Finite element matrix structure.
Definition: equations_matrices_routines.f90:72
base_routines
This module contains all the low-level base routines e.g., all debug, control, and low-level communic...
Definition: base_routines.f90:45
types::equations_matrices_linear_type
Contains information of the linear matrices for equations matrices.
Definition: types.f90:1479
iso_varying_string::trim
Definition: iso_varying_string.f90:204
equations_matrices_routines::equations_matrices_linear_storage_type_set
subroutine, public equations_matrices_linear_storage_type_set(EQUATIONS_MATRICES, STORAGE_TYPE, ERR, ERROR,)
Sets the storage type (sparsity) of the linear equations matrices.
Definition: equations_matrices_routines.f90:3898
equations_matrices_routines::equationsmatrices_linearstructuretypeset
subroutine, public equationsmatrices_linearstructuretypeset(EQUATIONS_MATRICES, STRUCTURE_TYPE, ERR, ERROR,)
Sets the structure (sparsity) of the linear equations matrices.
Definition: equations_matrices_routines.f90:4160
equations_mapping_routines::equations_mapping_create_finish
subroutine, public equations_mapping_create_finish(EQUATIONS_MAPPING, ERR, ERROR,)
Finishes the process of creating an equations mapping.
Definition: equations_mapping_routines.f90:982
control_loop_routines::control_loop_get
Returns the specified control loop as indexed by the control loop identifier from the control loop ro...
Definition: control_loop_routines.f90:86
equations_routines::equations_set_equations_get
subroutine, public equations_set_equations_get(EQUATIONS_SET, EQUATIONS, ERR, ERROR,)
Gets the equations for an equations set.
Definition: equations_routines.f90:897
types::equations_set_materials_type
Definition: types.f90:1881
constants::partial_derivative_first_derivative_map
integer(intg), dimension(4) partial_derivative_first_derivative_map
PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(nic) gives the partial derivative index for the first derivat...
Definition: constants.f90:254
equations_routines::equations_create_finish
subroutine, public equations_create_finish(EQUATIONS, ERR, ERROR,)
Finish the creation of equations.
Definition: equations_routines.f90:159
domain_mappings
This module handles all domain mappings routines.
Definition: domain_mappings.f90:45
problem_constants::problem_setup_finish_action
integer(intg), parameter problem_setup_finish_action
Finish setup action.
Definition: problem_constants.f90:274
maths::matrix_vector_product
Calculates and returns the matrix-vector-product A*b in the vector c.
Definition: maths.f90:203
equations_mapping_routines
This module handles all equations mapping routines.
Definition: equations_mapping_routines.f90:45
types::solver_equations_type
Contains information about the solver equations for a solver.
Definition: types.f90:2452
matrix_vector::matrix_compressed_row_storage_type
integer(intg), parameter, public matrix_compressed_row_storage_type
Matrix compressed row storage type.
Definition: matrix_vector.f90:167
equations_set_constants::equations_set_gfv_solution_method
integer(intg), parameter equations_set_gfv_solution_method
Grid-based Finite Volume solution method.
Definition: equations_set_constants.f90:354
equations_set_constants::equations_set_setup_geometry_type
integer(intg), parameter equations_set_setup_geometry_type
Geometry setup.
Definition: equations_set_constants.f90:323
constants::global_deriv_s1_s2
integer(intg), parameter global_deriv_s1_s2
Global Cross derivative in the s1 and s2 direction i.e., d^2u/ds1ds2.
Definition: constants.f90:216
types::problem_type
Contains information for a problem.
Definition: types.f90:3221
equations_set_constants::equations_set_classical_field_class
integer(intg), parameter equations_set_classical_field_class
Definition: equations_set_constants.f90:57
equations_set_constants::equations_set_hj_equation_type
integer(intg), parameter equations_set_hj_equation_type
Definition: equations_set_constants.f90:94
equations_set_constants::equations_linear
integer(intg), parameter equations_linear
The equations are linear.
Definition: equations_set_constants.f90:365
types::domain_nodes_type
Contains the topology information for the nodes of a domain.
Definition: types.f90:713
equations_matrices_routines::equations_matrices_create_finish
subroutine, public equations_matrices_create_finish(EQUATIONS_MATRICES, ERR, ERROR,)
Finishes the creation of the equations matrices and RHS for the the equations.
Definition: equations_matrices_routines.f90:404
distributed_matrix_vector
This module handles all distributed matrix vector routines.
Definition: distributed_matrix_vector.f90:45
constants::global_deriv_s1
integer(intg), parameter global_deriv_s1
First global derivative in the s1 direction i.e., du/ds1.
Definition: constants.f90:214
boundary_conditions_routines
This module handles all boundary conditions routines.
Definition: boundary_condition_routines.f90:45
solver_routines
This module handles all solver routines.
Definition: solver_routines.f90:45
equations_mapping_routines::equations_mapping_create_start
subroutine, public equations_mapping_create_start(EQUATIONS, EQUATIONS_MAPPING, ERR, ERROR,)
Finishes the process of creating an equations mapping for a equations set equations.
Definition: equations_mapping_routines.f90:1112
types::equations_matrix_type
Contains information about an equations matrix.
Definition: types.f90:1429
types::quadrature_scheme_type
Contains information for a particular quadrature scheme.
Definition: types.f90:141
matrix_vector
This module contains all routines dealing with (non-distributed) matrix and vectors types...
Definition: matrix_vector.f90:129
problem_constants::problem_generalised_hj_subtype
integer(intg), parameter problem_generalised_hj_subtype
Definition: problem_constants.f90:177
equations_routines::equations_linearity_type_set
subroutine, public equations_linearity_type_set(EQUATIONS, LINEARITY_TYPE, ERR, ERROR,)
Sets/changes the linearity type for equations.
Definition: equations_routines.f90:539
control_loop_routines::control_loop_create_start
subroutine, public control_loop_create_start(PROBLEM, CONTROL_LOOP, ERR, ERROR,)
Start the process of creating a control loop for a problem.
Definition: control_loop_routines.f90:194
problem_constants::problem_setup_solver_equations_type
integer(intg), parameter problem_setup_solver_equations_type
Solver equations setup for a problem.
Definition: problem_constants.f90:265
boundary_conditions_routines::boundary_conditions_set_local_dof
Sets a boundary condition on the specified local DOF.
Definition: boundary_condition_routines.f90:143
types::field_variable_type
Contains information for a field variable defined on a field.
Definition: types.f90:1289
equations_set_constants::equations_set_fd_solution_method
integer(intg), parameter equations_set_fd_solution_method
Finite Difference solution method.
Definition: equations_set_constants.f90:350
equations_matrices_routines::equations_matrices_sparse_matrices
integer(intg), parameter, public equations_matrices_sparse_matrices
Use sparse equations matrices.
Definition: equations_matrices_routines.f90:90
constants::global_deriv_s3
integer(intg), parameter global_deriv_s3
First global derivative in the s3 direction i.e., du/ds3.
Definition: constants.f90:217
types::equations_set_setup_type
Contains information on the setup information for an equations set.
Definition: types.f90:1866
types::domain_type
A pointer to the domain decomposition for this domain.
Definition: types.f90:938
iso_varying_string::index
Definition: iso_varying_string.f90:156
problem_constants::problem_setup_start_action
integer(intg), parameter problem_setup_start_action
Start setup action.
Definition: problem_constants.f90:273
solver_routines::solver_equations_time_dependence_type_set
subroutine, public solver_equations_time_dependence_type_set(SOLVER_EQUATIONS, TIME_DEPENDENCE_TYPE, ERR, ERROR,)
Sets/changes the time dependence type for solver equations.
Definition: solver_routines.f90:7788
control_loop_routines
This module handles all control loop routines.
Definition: control_loop_routines.f90:45
maths::matrix_product
Calculates and returns the matrix-product A*B in the matrix C.
Definition: maths.f90:167
boundary_conditions_routines::boundary_condition_fixed
integer(intg), parameter, public boundary_condition_fixed
The dof is fixed as a boundary condition.
Definition: boundary_condition_routines.f90:95
base_routines::errors
subroutine, public errors(NAME, ERR, ERROR)
Records the exiting error of the subroutine.
Definition: base_routines.f90:481
equations_set_constants
This module defines all constants shared across equations set routines.
Definition: equations_set_constants.f90:45
equations_set_constants::equations_set_bem_solution_method
integer(intg), parameter equations_set_bem_solution_method
Boundary Element Method solution method.
Definition: equations_set_constants.f90:349
solver_routines::solver_solver_equations_get
subroutine, public solver_solver_equations_get(SOLVER, SOLVER_EQUATIONS, ERR, ERROR,)
Returns a pointer to the solver equations for a solver.
Definition: solver_routines.f90:21230
problem_constants::problem_standard_hj_subtype
integer(intg), parameter problem_standard_hj_subtype
Definition: problem_constants.f90:176
equations_set_constants::equations_set_hj_equation_two_dim_2
integer(intg), parameter equations_set_hj_equation_two_dim_2
u=sin(x)sin(y)
Definition: equations_set_constants.f90:430
types::basis_type
Contains all information about a basis .
Definition: types.f90:184
equations_set_constants::equations_set_fv_solution_method
integer(intg), parameter equations_set_fv_solution_method
Finite Volume solution method.
Definition: equations_set_constants.f90:351
matrix_vector::matrix_block_storage_type
integer(intg), parameter, public matrix_block_storage_type
Matrix block storage type.
Definition: matrix_vector.f90:163
equations_set_constants::equations_set_setup_initial_type
integer(intg), parameter equations_set_setup_initial_type
Initial setup.
Definition: equations_set_constants.f90:322
control_loop_routines::control_loop_create_finish
recursive subroutine, public control_loop_create_finish(CONTROL_LOOP, ERR, ERROR,)
Finish the process of creating a control loop.
Definition: control_loop_routines.f90:153
constants::global_deriv_s1_s3
integer(intg), parameter global_deriv_s1_s3
Global Cross derivative in the s1 and s3 direction i.e., d^2u/ds1ds3.
Definition: constants.f90:218
equations_set_constants::equations_set_setup_analytic_type
integer(intg), parameter equations_set_setup_analytic_type
Analytic setup.
Definition: equations_set_constants.f90:329
base_routines::flagerror
Flags an error condition.
Definition: base_routines.f90:240
solver_routines::solver_linear_type
integer(intg), parameter, public solver_linear_type
A linear solver.
Definition: solver_routines.f90:88
types::equations_matrices_rhs_type
Contains information of the RHS vector for equations matrices.
Definition: types.f90:1500
constants::zero_tolerance
real(dp), parameter zero_tolerance
Definition: constants.f90:70
types::equations_mapping_linear_type
Contains information for mapping field variables to the linear matrices in the equations set of the m...
Definition: types.f90:1587
kinds
This module contains all kind definitions.
Definition: kinds.f90:45
equations_set_constants::equations_set_setup_finish_action
integer(intg), parameter equations_set_setup_finish_action
Finish setup action.
Definition: equations_set_constants.f90:339
input_output
This module handles all formating and input and output.
Definition: input_output.f90:45