OpenCMISS-Iron Internal API Documentation
basis_routines.f90
Go to the documentation of this file.
1 
43 
45 MODULE basis_routines
46 
47  USE base_routines
48  USE constants
49  USE input_output
50  USE iso_varying_string
51  USE kinds
52  USE strings
53  USE types
54 
55 #include "macros.h"
56 
57  IMPLICIT NONE
58 
59  PRIVATE
60 
61  !Module parameters
62 
68  INTEGER(INTG), PARAMETER :: basis_lagrange_hermite_tp_type=1
69  INTEGER(INTG), PARAMETER :: basis_simplex_type=2
70  INTEGER(INTG), PARAMETER :: basis_serendipity_type=3
71  INTEGER(INTG), PARAMETER :: basis_auxilliary_type=4
72  INTEGER(INTG), PARAMETER :: basis_b_spline_tp_type=5
73  INTEGER(INTG), PARAMETER :: basis_fourier_lagrange_hermite_tp_type=6
74  INTEGER(INTG), PARAMETER :: basis_extended_lagrange_tp_type=7
75  INTEGER(INTG), PARAMETER :: basis_radial_type=7
77 
82  INTEGER(INTG), PARAMETER :: basis_linear_lagrange_interpolation=1
83  INTEGER(INTG), PARAMETER :: basis_quadratic_lagrange_interpolation=2
84  INTEGER(INTG), PARAMETER :: basis_cubic_lagrange_interpolation=3
85  INTEGER(INTG), PARAMETER :: basis_cubic_hermite_interpolation=4
86  INTEGER(INTG), PARAMETER :: basis_quadratic1_hermite_interpolation=5
87  INTEGER(INTG), PARAMETER :: basis_quadratic2_hermite_interpolation=6
88  INTEGER(INTG), PARAMETER :: basis_linear_simplex_interpolation=7
89  INTEGER(INTG), PARAMETER :: basis_quadratic_simplex_interpolation=8
90  INTEGER(INTG), PARAMETER :: basis_cubic_simplex_interpolation=9
91  INTEGER(INTG), PARAMETER :: basis_gaussian_radial_interpolation=10
92  INTEGER(INTG), PARAMETER :: basis_multiquartic_radial_interpolation=11
94 
99  INTEGER(INTG), PARAMETER :: basis_lagrange_interpolation=1
100  INTEGER(INTG), PARAMETER :: basis_hermite_interpolation=2
101  INTEGER(INTG), PARAMETER :: basis_simplex_interpolation=3
102  INTEGER(INTG), PARAMETER :: basis_serendipity_interpolation=4
103  INTEGER(INTG), PARAMETER :: basis_transition_interpolation=5
104  INTEGER(INTG), PARAMETER :: basis_singular_interpolation=6
105  INTEGER(INTG), PARAMETER :: basis_fourier_interpolation=7
106  INTEGER(INTG), PARAMETER :: basis_radial_interpolation=8
108 
113  INTEGER(INTG), PARAMETER :: basis_linear_interpolation_order=1
114  INTEGER(INTG), PARAMETER :: basis_quadratic_interpolation_order=2
115  INTEGER(INTG), PARAMETER :: basis_cubic_interpolation_order=3
116  INTEGER(INTG), PARAMETER :: basis_quadratic1_interpolation_order=4
117  INTEGER(INTG), PARAMETER :: basis_quadratic2_interpolation_order=5
119 
124  INTEGER(INTG), PARAMETER :: basis_number_of_quadrature_scheme_types=4
125  INTEGER(INTG), PARAMETER :: basis_default_quadrature_scheme=1
126  INTEGER(INTG), PARAMETER :: basis_low_quadrature_scheme=2
127  INTEGER(INTG), PARAMETER :: basis_mid_quadrature_scheme=3
128  INTEGER(INTG), PARAMETER :: basis_high_quadrature_scheme=4
130 
135  INTEGER(INTG), PARAMETER :: basis_gauss_legendre_quadrature=1
136  INTEGER(INTG), PARAMETER :: basis_gauss_laguerre_quadrature=2
137  INTEGER(INTG), PARAMETER :: basis_guass_hermite_quadrature=3
138  INTEGER(INTG), PARAMETER :: basis_adaptive_gauss_legendre_quadrature=4
139  INTEGER(INTG), PARAMETER :: basis_gauss_simplex_quadrature=5
141 
146  INTEGER(INTG), PARAMETER :: basis_xi_collapsed=1
147  INTEGER(INTG), PARAMETER :: basis_collapsed_at_xi0=2
148  INTEGER(INTG), PARAMETER :: basis_collapsed_at_xi1=3
149  INTEGER(INTG), PARAMETER :: basis_not_collapsed=4
151 
152  !!Module types
153  !
154  !!>Contains information on the defined basis functions
155  !TYPE BASIS_FUNCTIONS_TYPE
156  ! PRIVATE
157  ! INTEGER(INTG) :: NUMBER_BASIS_FUNCTIONS !<The number of basis functions definegd
158  ! TYPE(BASIS_PTR_TYPE), POINTER :: BASES(:) !<The array of pointers to the defined basis functions
159  !END TYPE BASIS_FUNCTIONS_TYPE
160 
161  !Module variables
162 
163  TYPE(basis_functions_type) :: basis_functions
164 
165  !Interfaces
166 
168  INTERFACE basis_evaluate_xi
169  MODULE PROCEDURE basis_evaluate_xi_dp
170  END INTERFACE !BASIS_EVALUATE_XI
171 
173  INTERFACE basis_gauss_points_calculate
174  MODULE PROCEDURE basis_gauss_points_calculate_dp
175  END INTERFACE !BASIS_GAUSS_POINTS_CALCULATE
176 
178  INTERFACE basis_interpolate_gauss
179  MODULE PROCEDURE basis_interpolate_gauss_dp
180  END INTERFACE !BASIS_INTERPOLATE_GAUSS
181 
183  INTERFACE basis_interpolate_xi
184  MODULE PROCEDURE basis_interpolate_xi_dp
185  END INTERFACE !BASIS_INTERPOLATE_XI
186 
188  INTERFACE basis_interpolate_local_face_gauss
189  MODULE PROCEDURE basis_interpolate_local_face_gauss_dp
190  END INTERFACE !BASIS_INTERPOLATE_LOCAL_FACE_GAUSS
191 
193  INTERFACE basis_interpolation_xi_set
194  MODULE PROCEDURE basis_interpolation_xi_set_number
195  MODULE PROCEDURE basis_interpolation_xi_set_ptr
196  END INTERFACE !BASIS_INTERPOLATION_XI_SET
197 
199  INTERFACE basis_number_of_xi_set
200  MODULE PROCEDURE basis_number_of_xi_set_number
201  MODULE PROCEDURE basis_number_of_xi_set_ptr
202  END INTERFACE !BASIS_NUMBER_OF_XI_SET
203 
205  INTERFACE basis_type_set
206  MODULE PROCEDURE basis_type_set_number
207  MODULE PROCEDURE basis_type_set_ptr
208  END INTERFACE !BASIS_TYPE_SET
209 
211  INTERFACE basis_collapsed_xi_set
212  MODULE PROCEDURE basis_collapsed_xi_set_number
213  MODULE PROCEDURE basis_collapsed_xi_set_ptr
214  END INTERFACE !BASIS_COLLAPSED_XI_SET
215 
217  INTERFACE basis_quadrature_order_set
218  MODULE PROCEDURE basis_quadrature_order_set_number
219  MODULE PROCEDURE basis_quadrature_order_set_ptr
220  END INTERFACE !BASIS_QUADRATURE_ORDER_SET
221 
223  INTERFACE basis_quadrature_type_set
224  MODULE PROCEDURE basis_quadrature_type_set_number
225  MODULE PROCEDURE basis_quadrature_type_set_ptr
226  END INTERFACE !BASIS_QUADRATURE_TYPE_SET
227 
228  INTERFACE simplex_linear_evaluate
229  MODULE PROCEDURE simplex_linear_evaluate_dp
230  END INTERFACE !SIMPLEX_LINEAR_EVALUATE
231 
232  INTERFACE simplex_quadratic_evaluate
233  MODULE PROCEDURE simplex_quadratic_evaluate_dp
234  END INTERFACE !SIMPLEX_QUADRATIC_EVALUATE
235 
236  INTERFACE simplex_cubic_evaluate
237  MODULE PROCEDURE simplex_cubic_evaluate_dp
238  END INTERFACE !SIMPLEX_CUBIC_EVALUATE
239 
241  INTERFACE basis_lhtp_basis_evaluate
242  MODULE PROCEDURE basis_lhtp_basis_evaluate_dp
243  END INTERFACE !BASIS_LHTP_BASIS_EVALUATE
244 
245  PUBLIC basis_lagrange_hermite_tp_type,basis_simplex_type,basis_serendipity_type,basis_auxilliary_type,basis_b_spline_tp_type, &
246  & basis_fourier_lagrange_hermite_tp_type,basis_extended_lagrange_tp_type, basis_radial_type
247 
248  PUBLIC basis_linear_lagrange_interpolation,basis_quadratic_lagrange_interpolation,basis_cubic_lagrange_interpolation, &
249  & basis_cubic_hermite_interpolation,basis_quadratic1_hermite_interpolation,basis_quadratic2_hermite_interpolation, &
250  & basis_linear_simplex_interpolation,basis_quadratic_simplex_interpolation,basis_cubic_simplex_interpolation, &
251  & basis_gaussian_radial_interpolation, basis_multiquartic_radial_interpolation
252 
253  PUBLIC basis_linear_interpolation_order,basis_quadratic_interpolation_order,basis_cubic_interpolation_order, &
254  & basis_quadratic1_interpolation_order,basis_quadratic2_interpolation_order
255 
256  PUBLIC basis_default_quadrature_scheme,basis_low_quadrature_scheme,basis_mid_quadrature_scheme,basis_high_quadrature_scheme
257 
258  PUBLIC basis_gauss_legendre_quadrature,basis_gauss_laguerre_quadrature,basis_guass_hermite_quadrature,&
259  & basis_adaptive_gauss_legendre_quadrature,basis_gauss_simplex_quadrature
260 
261  PUBLIC basis_xi_collapsed,basis_collapsed_at_xi0,basis_collapsed_at_xi1,basis_not_collapsed, basis_functions
262 
263  PUBLIC basis_collapsed_xi_set
264 
265  PUBLIC basis_evaluate_xi
266 
267  PUBLIC basis_gauss_points_calculate
268 
269  PUBLIC basis_interpolate_gauss,basis_interpolate_xi,basis_interpolate_local_face_gauss
270 
271  PUBLIC basis_local_node_xi_calculate
272 
273  PUBLIC basis_number_of_local_nodes_get
274 
275  PUBLIC basis_interpolation_xi_set
276 
277  PUBLIC basis_number_of_xi_set
278 
279  PUBLIC basis_quadrature_number_of_gauss_xi_set
280 
281  PUBLIC basis_quadrature_destroy,basis_quadrature_order_set,basis_quadrature_type_set
282 
283  PUBLIC basis_type_set,basis_quadraturelocalfacegaussevaluateset
284 
285  PUBLIC basis_create_start,basis_create_finish
286 
287  PUBLIC basis_destroy
288 
289  PUBLIC bases_finalise,bases_initialise
290 
291  PUBLIC basis_user_number_find
292 
293  PUBLIC basis_collapsed_xi_get
294 
295  PUBLIC basis_interpolation_xi_get
296 
297  PUBLIC basis_number_of_xi_get
298 
299  PUBLIC basis_quadrature_number_of_gauss_xi_get
300 
301  PUBLIC basis_quadrature_single_gauss_xi_get, basis_quadrature_multiple_gauss_xi_get
302 
303  PUBLIC basis_quadrature_order_get
304 
305  PUBLIC basis_quadrature_type_get
306 
307  PUBLIC basis_type_get
308 
309 
310 CONTAINS
311 
312  !
313  !================================================================================================================================
314  !
315 
317  SUBROUTINE bases_finalise(ERR,ERROR,*)
318 
319  !Argument variables
320  INTEGER(INTG), INTENT(OUT) :: ERR
321  TYPE(varying_string), INTENT(OUT) :: ERROR
322  !Local Variables
323 
324  enters("BASES_FINALISE",err,error,*999)
325 
326  !Destroy any created basis functions
327  DO WHILE(basis_functions%NUMBER_BASIS_FUNCTIONS>0)
328  CALL basis_destroy(basis_functions%BASES(1)%PTR,err,error,*999)
329  ENDDO !nb
330  !Destroy basis functions and deallocated any memory allocated
331  basis_functions%NUMBER_BASIS_FUNCTIONS=0
332  IF(ASSOCIATED(basis_functions%BASES)) DEALLOCATE(basis_functions%BASES)
333 
334  exits("BASES_FINALISE")
335  RETURN
336 999 errorsexits("BASES_FINALISE",err,error)
337  RETURN 1
338 
339  END SUBROUTINE bases_finalise
340 
341  !
342  !================================================================================================================================
343  !
344 
346  SUBROUTINE bases_initialise(ERR,ERROR,*)
347 
348  !Argument variables
349  INTEGER(INTG), INTENT(OUT) :: ERR
350  TYPE(varying_string), INTENT(OUT) :: ERROR
351  !Local Variables
352 
353  enters("BASES_INITIALISE",err,error,*999)
354 
355  basis_functions%NUMBER_BASIS_FUNCTIONS=0
356  NULLIFY(basis_functions%BASES)
357 
358  exits("BASES_INITIALISE")
359  RETURN
360 999 errorsexits("BASES_INITIALISE",err,error)
361  RETURN 1
362 
363  END SUBROUTINE bases_initialise
364 
365  !
366  !================================================================================================================================
367  !
368 
370  SUBROUTINE basis_create_finish(BASIS,ERR,ERROR,*)
371 
372  !Argument variables
373  TYPE(basis_type), POINTER :: BASIS
374  INTEGER(INTG), INTENT(OUT) :: ERR
375  TYPE(varying_string), INTENT(OUT) :: ERROR
376  !Local Variables
377  INTEGER(INTG) :: ni,nic,nn,nn1,nn2,nn3,nn4,ns,local_line_idx,local_face_idx
378  TYPE(varying_string) :: LOCAL_ERROR
379 
380  enters("BASIS_CREATE_FINISH",err,error,*999)
381 
382  IF(ASSOCIATED(basis)) THEN
383  SELECT CASE(basis%TYPE)
384  CASE(basis_lagrange_hermite_tp_type)
385  CALL basis_lhtp_family_create(basis,err,error,*999)
386  CASE(basis_simplex_type)
387  CALL basis_simplex_family_create(basis,err,error,*999)
388  CASE(basis_radial_type)
389  CALL basis_radial_family_create(basis,err,error,*999)
390  CASE DEFAULT
391  local_error="Basis type "//trim(number_to_vstring(basis%TYPE,"*",err,error))//" is invalid or not implemented"
392  CALL flagerror(local_error,err,error,*999)
393  END SELECT
394  basis%BASIS_FINISHED=.true.
395  ELSE
396  CALL flagerror("Basis is not associated",err,error,*999)
397  ENDIF
398 
399  IF(diagnostics1) THEN
400  CALL write_string_value(diagnostic_output_type,"Basis user number = ",basis%USER_NUMBER,err,error,*999)
401  CALL write_string_value(diagnostic_output_type," Basis family number = ",basis%FAMILY_NUMBER,err,error,*999)
402  CALL write_string_value(diagnostic_output_type," Basis global number = ",basis%GLOBAL_NUMBER,err,error,*999)
403  CALL write_string_value(diagnostic_output_type," Basis type = ",basis%TYPE,err,error,*999)
404  CALL write_string_value(diagnostic_output_type," Basis degenerate = ",basis%DEGENERATE,err,error,*999)
405  CALL write_string_value(diagnostic_output_type," Number of Xi directions = ",basis%NUMBER_OF_XI,err,error,*999)
406  CALL write_string_value(diagnostic_output_type," Number of Xi coordinates = ",basis%NUMBER_OF_XI_COORDINATES,err,error,*999)
407 
408  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_XI_COORDINATES,4,4,basis%INTERPOLATION_TYPE, &
409  & '(" Interpolation type(nic):",4(X,I2))','(25X,4(X,I2))',err,error,*999)
410  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_XI_COORDINATES,4,4,basis%INTERPOLATION_ORDER, &
411  & '(" Interpolation order(nic):",4(X,I2))','(26X,4(X,I2))',err,error,*999)
412  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_XI,3,3,basis%COLLAPSED_XI, &
413  & '(" Collapsed xi(ni):",3(X,I2))','(26X,3(X,I2))',err,error,*999)
414  CALL write_string_value(diagnostic_output_type," Number of partial derivatives = ",basis%NUMBER_OF_PARTIAL_DERIVATIVES, &
415  & err,error,*999)
416  CALL write_string_value(diagnostic_output_type," Total number of nodes = ",basis%NUMBER_OF_NODES,err,error,*999)
417  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_XI_COORDINATES,4,4,basis%NUMBER_OF_NODES_XIC, &
418  & '(" Number of nodes(nic):",4(X,I2))','(22X,4(X,I2))',err,error,*999)
419  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_NODES,8,8,basis%NUMBER_OF_DERIVATIVES, &
420  & '(" Number of derivatives(nn):",8(X,I2))','(28X,8(X,I2))',err,error,*999)
421  CALL write_string_value(diagnostic_output_type," Number of element parameters = ",basis%NUMBER_OF_ELEMENT_PARAMETERS, &
422  & err,error,*999)
423 ! CPB 23/07/07 Doxygen may or may not like this line for some reason????
424  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_NODES,8,8,basis%NODE_AT_COLLAPSE, &
425  & '(" Node at collapse(nn):",8(X,L1))','(23X,8(X,L1))',err,error,*999)
426  CALL write_string(diagnostic_output_type," Node position index:",err,error,*999)
427  DO nic=1,basis%NUMBER_OF_XI_COORDINATES
428  CALL write_string_value(diagnostic_output_type," Xic = ",nic,err,error,*999)
429  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_NODES,16,16,basis%NODE_POSITION_INDEX(:,nic), &
430  & '(" INDEX(nn) :",16(X,I2))','(18X,16(X,I2))',err,error,*999)
431  ENDDO !ni
432 
433  CALL write_string(diagnostic_output_type," Inverse node position index:",err,error,*999)
434  SELECT CASE(basis%NUMBER_OF_XI_COORDINATES)
435  CASE(1)
436  DO nn1=1,basis%NUMBER_OF_NODES_XIC(1)
437  CALL write_string_value(diagnostic_output_type," Xic = 1, Node position index = ",nn1,err,error,*999)
438  CALL write_string_value(diagnostic_output_type," INDEX = ",basis%NODE_POSITION_INDEX_INV(nn1,1,1,1), &
439  & err,error,*999)
440  ENDDO !nn1
441  CASE(2)
442  DO nn2=1,basis%NUMBER_OF_NODES_XIC(2)
443  CALL write_string_value(diagnostic_output_type," Xic = 2, Node position index = ",nn2,err,error,*999)
444  DO nn1=1,basis%NUMBER_OF_NODES_XIC(1)
445  CALL write_string_value(diagnostic_output_type," Xic = 1, Node position index = ",nn1,err,error,*999)
446  CALL write_string_value(diagnostic_output_type," INDEX = ",basis%NODE_POSITION_INDEX_INV(nn1,nn2,1,1), &
447  & err,error,*999)
448  ENDDO !nn1
449  ENDDO !nn2
450  CASE(3)
451  DO nn3=1,basis%NUMBER_OF_NODES_XIC(3)
452  CALL write_string_value(diagnostic_output_type," Xic = 3, Node position index = ",nn3,err,error,*999)
453  DO nn2=1,basis%NUMBER_OF_NODES_XIC(2)
454  CALL write_string_value(diagnostic_output_type," Xic = 2, Node position index = ",nn2,err,error,*999)
455  DO nn1=1,basis%NUMBER_OF_NODES_XIC(1)
456  CALL write_string_value(diagnostic_output_type," Xic = 1, Node position index = ",nn1,err,error,*999)
457  CALL write_string_value(diagnostic_output_type," INDEX = ",basis%NODE_POSITION_INDEX_INV(nn1,nn2,nn3,1), &
458  & err,error,*999)
459  ENDDO !nn1
460  ENDDO !nn2
461  ENDDO !nn3
462  CASE(4)
463  DO nn4=1,basis%NUMBER_OF_NODES_XIC(4)
464  CALL write_string_value(diagnostic_output_type," Xic = 4, Node position index = ",nn4,err,error,*999)
465  DO nn3=1,basis%NUMBER_OF_NODES_XIC(3)
466  CALL write_string_value(diagnostic_output_type," Xic = 3, Node position index = ",nn3,err,error,*999)
467  DO nn2=1,basis%NUMBER_OF_NODES_XIC(2)
468  CALL write_string_value(diagnostic_output_type," Xic = 2, Node position index = ",nn2,err,error,*999)
469  DO nn1=1,basis%NUMBER_OF_NODES_XIC(1)
470  CALL write_string_value(diagnostic_output_type," Xic = 1, Node position index = ",nn1,err,error,*999)
471  CALL write_string_value(diagnostic_output_type," INDEX = " &
472  & ,basis%NODE_POSITION_INDEX_INV(nn1,nn2,nn3,nn4),err,error,*999)
473  ENDDO !nn1
474  ENDDO !nn2
475  ENDDO !nn3
476  ENDDO !nn4
477  CASE DEFAULT
478  CALL flagerror("Invalid number of xi coordinates",err,error,*999)
479  END SELECT
480  CALL write_string(diagnostic_output_type," Derivative order index:",err,error,*999)
481  DO ni=1,basis%NUMBER_OF_XI
482  CALL write_string_value(diagnostic_output_type," Xi = ",ni,err,error,*999)
483  DO nn=1,basis%NUMBER_OF_NODES
484  CALL write_string_value(diagnostic_output_type," Node = ",nn,err,error,*999)
485  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_DERIVATIVES(nn),8,8, &
486  & basis%DERIVATIVE_ORDER_INDEX(:,nn,ni),'(" INDEX(nk):",8(X,I2))','(18X,8(X,I2))',err,error,*999)
487  ENDDO !nn
488  ENDDO !ni
489  CALL write_string(diagnostic_output_type," Inverse derivative order index:",err,error,*999)
490  CALL write_string(diagnostic_output_type," Element parameter index:",err,error,*999)
491  DO nn=1,basis%NUMBER_OF_NODES
492  CALL write_string_value(diagnostic_output_type," Node = ",nn,err,error,*999)
493  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_DERIVATIVES(nn),8,8, &
494  & basis%ELEMENT_PARAMETER_INDEX(:,nn),'(" INDEX(nk) :",8(X,I2))','(18X,8(X,I2))',err,error,*999)
495  ENDDO !nn
496  CALL write_string(diagnostic_output_type," Inverse element parameter index:",err,error,*999)
497  DO ns=1,basis%NUMBER_OF_ELEMENT_PARAMETERS
498  CALL write_string_value(diagnostic_output_type," Element parameter = ",ns,err,error,*999)
499  CALL write_string_vector(diagnostic_output_type,1,1,2,2,2, &
500  & basis%ELEMENT_PARAMETER_INDEX_INV(:,ns),'(" INDEX(:) :",2(X,I2))','(18X,2(X,I2))',err,error,*999)
501  ENDDO !ns
502  CALL write_string(diagnostic_output_type," Partial derivative index:",err,error,*999)
503  DO nn=1,basis%NUMBER_OF_NODES
504  CALL write_string_value(diagnostic_output_type," Node = ",nn,err,error,*999)
505  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_DERIVATIVES(nn),8,8, &
506  & basis%PARTIAL_DERIVATIVE_INDEX(:,nn),'(" INDEX(nk) :",8(X,I2))','(18X,8(X,I2))',err,error,*999)
507  ENDDO !nn
508  IF(basis%NUMBER_OF_XI==3) THEN
509  CALL write_string(diagnostic_output_type," Local faces:",err,error,*999)
510  CALL write_string_value(diagnostic_output_type," Number of local faces = ",basis%NUMBER_OF_LOCAL_FACES,err,error,*999)
511  DO local_face_idx=1,basis%NUMBER_OF_LOCAL_FACES
512  CALL write_string_value(diagnostic_output_type," Local face = ",local_face_idx,err,error,*999)
513  CALL write_string_value(diagnostic_output_type," Local face xi direction = ", &
514  & basis%LOCAL_FACE_XI_DIRECTION(local_face_idx),err,error,*999)
515  CALL write_string_value(diagnostic_output_type," Number of nodes in local face = ", &
516  & basis%NUMBER_OF_NODES_IN_LOCAL_FACE(local_face_idx),err,error,*999)
517  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_NODES_IN_LOCAL_FACE(local_face_idx),4,4, &
518  & basis%NODE_NUMBERS_IN_LOCAL_FACE(:,local_face_idx),'(" Nodes in local face :",4(X,I2))','(33X,4(X,I2))', &
519  & err,error,*999)
520  DO nn=1,basis%NUMBER_OF_NODES_IN_LOCAL_FACE(local_face_idx)
521  CALL write_string_value(diagnostic_output_type," Local face node: ",nn,err,error,*999)
522  CALL write_string_vector(diagnostic_output_type,1,1,basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0,nn,local_face_idx),4,4, &
523  & basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(1:,nn,local_face_idx),'(" Derivatives in local face :",4(X,I2))', &
524  & '(33X,4(X,I2))',err,error,*999)
525  ENDDO
526  ENDDO !ni
527  ENDIF
528  CALL write_string(diagnostic_output_type," Local lines:",err,error,*999)
529  CALL write_string_value(diagnostic_output_type," Number of local lines = ",basis%NUMBER_OF_LOCAL_LINES,err,error,*999)
530  DO local_line_idx=1,basis%NUMBER_OF_LOCAL_LINES
531  CALL write_string_value(diagnostic_output_type," Local line = ",local_line_idx,err,error,*999)
532  CALL write_string_value(diagnostic_output_type," Local line xi direction = ", &
533  & basis%LOCAL_LINE_XI_DIRECTION(local_line_idx),err,error,*999)
534  CALL write_string_value(diagnostic_output_type," Number of nodes in local line = ", &
535  & basis%NUMBER_OF_NODES_IN_LOCAL_LINE(local_line_idx),err,error,*999)
536  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_NODES_IN_LOCAL_LINE(local_line_idx),4,4, &
537  & basis%NODE_NUMBERS_IN_LOCAL_LINE(:,local_line_idx),'(" Nodes in local line :",4(X,I2))','(33X,4(X,I2))', &
538  & err,error,*999)
539  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_NODES_IN_LOCAL_LINE(local_line_idx),4,4, &
540  & basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(:,local_line_idx),'(" Derivatives in local line :",4(X,I2))', &
541  & '(33X,4(X,I2))',err,error,*999)
542  IF(basis%NUMBER_OF_XI==2) THEN
543  CALL write_string_value(diagnostic_output_type," Local line xi normal = ", &
544  & basis%LOCAL_XI_NORMAL(local_line_idx),err,error,*999)
545  ENDIF
546  ENDDO !ni
547  CALL write_string_value(diagnostic_output_type," Number of sub-bases = ",basis%NUMBER_OF_SUB_BASES,err,error,*999)
548  ENDIF
549 
550  exits("BASIS_CREATE_FINISH")
551  RETURN
552 999 errorsexits("BASIS_CREATE_FINISH",err,error)
553  RETURN 1
554 
555  END SUBROUTINE basis_create_finish
556 
557  !
558  !================================================================================================================================
559  !
560 
574  SUBROUTINE basis_create_start(USER_NUMBER,BASIS,ERR,ERROR,*)
575 
576  !Argument variables
577  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
578  TYPE(basis_type), POINTER :: BASIS
579  INTEGER(INTG), INTENT(OUT) :: ERR
580  TYPE(varying_string), INTENT(OUT) :: ERROR
581  !Local Variables
582  INTEGER(INTG) :: nb
583  TYPE(basis_type), POINTER :: NEW_BASIS
584  TYPE(basis_ptr_type), POINTER :: NEW_BASES(:)
585 
586  NULLIFY(new_basis)
587  NULLIFY(new_bases)
588 
589  enters("BASIS_CREATE_START",err,error,*999)
590 
591  IF(ASSOCIATED(basis)) THEN
592  CALL flagerror("Basis is already associated",err,error,*999)
593  ELSE
594  !See if basis number has already been created
595  CALL basis_user_number_find(user_number,basis,err,error,*999)
596  IF(ASSOCIATED(basis)) THEN
597  CALL flagerror("Basis number is already defined",err,error,*999)
598  ELSE
599  !Allocate new basis function and add it to the basis functions
600  ALLOCATE(new_bases(basis_functions%NUMBER_BASIS_FUNCTIONS+1),stat=err)
601  IF(err/=0) CALL flagerror("Could not allocate NEW_BASES",err,error,*999)
602  ALLOCATE(new_basis,stat=err)
603  IF(err/=0) CALL flagerror("Could not allocate NEW_BASIS",err,error,*999)
604  CALL basis_initialise(new_basis,err,error,*999)
605  DO nb=1,basis_functions%NUMBER_BASIS_FUNCTIONS
606  new_bases(nb)%PTR=>basis_functions%BASES(nb)%PTR
607  ENDDO !nb
608  basis_functions%NUMBER_BASIS_FUNCTIONS=basis_functions%NUMBER_BASIS_FUNCTIONS+1
609  new_bases(basis_functions%NUMBER_BASIS_FUNCTIONS)%PTR=>new_basis
610  IF(ASSOCIATED(basis_functions%BASES)) DEALLOCATE(basis_functions%BASES)
611  basis_functions%BASES=>new_bases
612  !Initialise the new basis function pointers
613  new_basis%NUMBER_OF_SUB_BASES=0
614  NULLIFY(new_basis%SUB_BASES)
615  NULLIFY(new_basis%PARENT_BASIS)
616  !Set the basis parameters
617  new_basis%BASIS_FINISHED=.false.
618  new_basis%USER_NUMBER=user_number
619  new_basis%FAMILY_NUMBER=0
620  new_basis%GLOBAL_NUMBER=basis_functions%NUMBER_BASIS_FUNCTIONS
621  !Set the default basis parameters
622  new_basis%TYPE=basis_lagrange_hermite_tp_type
623  new_basis%NUMBER_OF_XI=3
624  ALLOCATE(new_basis%INTERPOLATION_XI(3),stat=err)
625  IF(err/=0) CALL flagerror("Could not allocate basis interpolation xi",err,error,*999)
626  new_basis%INTERPOLATION_XI(1:3)=[basis_linear_lagrange_interpolation,basis_linear_lagrange_interpolation, &
627  & basis_linear_lagrange_interpolation]
628  ALLOCATE(new_basis%COLLAPSED_XI(3),stat=err)
629  IF(err/=0) CALL flagerror("Could not allocate basis collapsed xi",err,error,*999)
630  new_basis%COLLAPSED_XI(1:3)=basis_not_collapsed
631  !Initialise the basis quadrature
632  NULLIFY(new_basis%QUADRATURE%BASIS)
633  CALL basis_quadrature_initialise(new_basis,err,error,*999)
634 
635  basis=>new_basis
636  ENDIF
637  ENDIF
638 
639  exits("BASIS_CREATE_START")
640  RETURN
641 999 IF(ASSOCIATED(new_basis)) CALL basis_destroy(new_basis,err,error,*998)
642 998 IF(ASSOCIATED(new_bases)) DEALLOCATE(new_bases)
643  NULLIFY(basis)
644  errorsexits("BASIS_CREATE_START",err,error)
645  RETURN 1
646 
647  END SUBROUTINE basis_create_start
648 
649  !
650  !================================================================================================================================
651  !
652 
654  RECURSIVE SUBROUTINE basis_destroy_number(USER_NUMBER,ERR,ERROR,*)
655 
656  !Argument variables
657  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
658  INTEGER(INTG), INTENT(OUT) :: ERR
659  TYPE(varying_string), INTENT(OUT) :: ERROR
660  !Local Variables
661 
662  enters("BASIS_DESTROY_NUMBER",err,error,*999)
663 
664  CALL basis_family_destroy(user_number,0,err,error,*999)
665 
666  exits("BASIS_DESTROY_NUMBER")
667  RETURN
668 999 errorsexits("BASIS_DESTROY_NUMBER",err,error)
669  RETURN 1
670 
671  END SUBROUTINE basis_destroy_number
672 
673  !
674  !================================================================================================================================
675  !
676 
678  RECURSIVE SUBROUTINE basis_destroy(BASIS,ERR,ERROR,*)
679 
680  !Argument variables
681  TYPE(basis_type), POINTER :: BASIS
682  INTEGER(INTG), INTENT(OUT) :: ERR
683  TYPE(varying_string), INTENT(OUT) :: ERROR
684  !Local Variables
685  INTEGER(INTG) :: USER_NUMBER
686 
687  enters("BASIS_DESTROY",err,error,*999)
688 
689  IF(ASSOCIATED(basis)) THEN
690  user_number=basis%USER_NUMBER
691  CALL basis_family_destroy(user_number,0,err,error,*999)
692  !NULLIFY(BASIS)
693  ELSE
694  CALL flagerror("Basis is not associated.",err,error,*999)
695  ENDIF
696 
697  exits("BASIS_DESTROY")
698  RETURN
699 999 errorsexits("BASIS_DESTROY",err,error)
700  RETURN 1
701 
702  END SUBROUTINE basis_destroy
703 
704  !
705  !================================================================================================================================
706  !
707 
710  FUNCTION basis_evaluate_xi_dp(BASIS,ELEMENT_PARAMETER_INDEX,PARTIAL_DERIV_INDEX,XI,ERR,ERROR)
711 
712  !Argument variables
713  TYPE(basis_type), POINTER :: BASIS
714  INTEGER(INTG), INTENT(IN) :: ELEMENT_PARAMETER_INDEX
715  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIV_INDEX
716  REAL(DP), INTENT(IN) :: XI(:)
717  INTEGER(INTG), INTENT(OUT) :: ERR
718  TYPE(varying_string), INTENT(OUT) :: ERROR
719  !Function variable
720  REAL(DP) :: BASIS_EVALUATE_XI_DP
721  !Local Variables
722  INTEGER(INTG) :: nn,nk
723  REAL(DP) :: XIL(size(xi,1)+1)
724  TYPE(varying_string) :: LOCAL_ERROR
725 
726  enters("BASIS_EVALUATE_XI_DP",err,error,*999)
727 
728  basis_evaluate_xi_dp=0.0_dp
729  IF(ASSOCIATED(basis)) THEN
730  IF(element_parameter_index>0.AND.element_parameter_index<=basis%NUMBER_OF_ELEMENT_PARAMETERS) THEN
731  SELECT CASE(basis%TYPE)
732  CASE(basis_lagrange_hermite_tp_type)
733  nn=basis%ELEMENT_PARAMETER_INDEX_INV(1,element_parameter_index)
734  nk=basis%ELEMENT_PARAMETER_INDEX_INV(2,element_parameter_index)
735  basis_evaluate_xi_dp=basis_lhtp_basis_evaluate(basis,nn,nk,partial_deriv_index,xi,err,error)
736  IF(err/=0) GOTO 999
737  CASE(basis_simplex_type)
738  !Create the area coordinates from the xi coordinates
739  xil(1:SIZE(xi,1))=1.0_dp-xi
740  xil(SIZE(xi,1)+1)=sum(xi)-(SIZE(xi,1)-1.0_dp)
741  nn=basis%ELEMENT_PARAMETER_INDEX_INV(1,element_parameter_index)
742  basis_evaluate_xi_dp=basis_simplex_basis_evaluate(basis,nn,partial_deriv_index,xil,err,error)
743  IF(err/=0) GOTO 999
744  CASE(basis_radial_type)
745 
746 
747  IF(err/=0) GOTO 999
748  CASE DEFAULT
749  local_error="Basis type "//trim(number_to_vstring(basis%TYPE,"*",err,error))//" is invalid or not implemented."
750  CALL flagerror(local_error,err,error,*999)
751  END SELECT
752  ELSE
753  local_error="The specified element parameter index of "// &
754  & trim(number_to_vstring(element_parameter_index,"*",err,error))// &
755  & " is invalid. The index must be > 0 and <= "// &
756  & trim(number_to_vstring(basis%NUMBER_OF_ELEMENT_PARAMETERS,"*",err,error))//"."
757  CALL flagerror(local_error,err,error,*999)
758  ENDIF
759  ELSE
760  CALL flagerror("Basis is not associated.",err,error,*999)
761  ENDIF
762 
763  exits("BASIS_EVALUATE_XI_DP")
764  RETURN
765 999 errorsexits("BASIS_EVALUATE_XI_DP",err,error)
766  RETURN
767 
768  END FUNCTION basis_evaluate_xi_dp
769 
770  !
771  !================================================================================================================================
772  !
773 
775  SUBROUTINE basis_gauss_points_calculate_dp(basis,order,numCoords,numberGaussPoints,gaussPoints,gaussWeights,err,error,*)
776 
777  !Argument variables
778  TYPE(basis_type), POINTER :: basis
779  INTEGER(INTG), INTENT(IN) :: order
780  INTEGER(INTG), INTENT(IN) :: numCoords
781  INTEGER(INTG), INTENT(OUT) :: numberGaussPoints
782  REAL(DP), ALLOCATABLE, INTENT(OUT) :: gaussPoints(:,:)
783  REAL(DP), ALLOCATABLE, INTENT(OUT) :: gaussWeights(:)
784  INTEGER(INTG), INTENT(OUT) :: err
785  TYPE(varying_string), INTENT(OUT) :: error
786  !Local Variables
787  INTEGER(INTG) :: number_of_vertices,nj,ng,i,j,k,NUM_GAUSS_1,NUM_GAUSS_2,NUM_GAUSS_3,MAX_GAUSS
788  REAL(DP), ALLOCATABLE :: XICOORDS(:,:),W(:,:),XI_MATRIX(:,:,:,:),XI(:)
789  TYPE(varying_string) :: LOCAL_ERROR
790 
791  enters("BASIS_GAUSS_POINTS_CALCULATE_DP",err,error,*999)
792  number_of_vertices=0
793  IF(ASSOCIATED(basis)) THEN
794  !current code assumes same order in each direction
795  SELECT CASE(numcoords)
796  CASE(1)
797  num_gauss_1=order
798  num_gauss_2=1
799  num_gauss_3=1
800  max_gauss=order
801  CASE(2)
802  num_gauss_1=order
803  num_gauss_2=order
804  num_gauss_3=1
805  max_gauss=order*order
806  CASE(3)
807  num_gauss_1=order
808  num_gauss_2=order
809  num_gauss_3=order
810  max_gauss=order*order*order
811  CASE DEFAULT
812  local_error="Number of coordinates " &
813  & //trim(number_to_vstring(numcoords,"*",err,error))//" is invalid or not implemented"
814  CALL flagerror(local_error,err,error,*999)
815  END SELECT
816  !Allocate arrays
817  ALLOCATE(w(max_gauss,numcoords),stat=err)
818  IF(err/=0) CALL flagerror("Could not allocate weights",err,error,*999)
819  ALLOCATE(xi(numcoords),stat=err)
820  IF(err/=0) CALL flagerror("Could not allocate gauss point coordinates",err,error,*999)
821  ALLOCATE(xicoords(max_gauss,numcoords),stat=err)
822  IF(err/=0) CALL flagerror("Could not allocate gauss point coordinates",err,error,*999)
823  ALLOCATE(xi_matrix(max_gauss,max_gauss,max_gauss,numcoords),stat=err)
824  IF(err/=0) CALL flagerror("Could not allocate XI matrix",err,error,*999)
825 
826  SELECT CASE(basis%TYPE)
827  CASE(basis_lagrange_hermite_tp_type)
828  ALLOCATE(gausspoints(numcoords,max_gauss),stat=err)
829  IF(err/=0) CALL flagerror("Could not allocate weights",err,error,*999)
830  ALLOCATE(gaussweights(max_gauss),stat=err)
831  IF(err/=0) CALL flagerror("Could not allocate weights",err,error,*999)
832  DO nj=1,numcoords
833  CALL gauss_legendre(order,0.0_dp,1.0_dp,xicoords(1:order,nj),w(1:order,nj),err,error,*999)
834  IF(err/=0) GOTO 999
835  ENDDO
836  !Form gauss point array for lagrange hermite type.
837  numbergausspoints=0
838  DO k=1,num_gauss_3
839  DO j=1,num_gauss_2
840  DO i=1,num_gauss_1
841  xi_matrix(i,j,k,1)=xicoords(i,1)
842  xi_matrix(i,j,k,2)=xicoords(j,2)
843  xi_matrix(i,j,k,3)=xicoords(k,3)
844  xi(1:numcoords)=xi_matrix(i,j,k,1:numcoords)
845  ng=i+(j-1+(k-1)*num_gauss_2)*num_gauss_1
846  gaussweights(ng)=w(i,1)*w(j,2)*w(k,3)
847  gausspoints(1:numcoords,ng)=xi(1:numcoords)
848  numbergausspoints=numbergausspoints+1
849  ENDDO
850  ENDDO
851  ENDDO
852  CASE(basis_simplex_type)
853  IF(numcoords==1) THEN
854  number_of_vertices=2
855  ELSEIF(numcoords==2) THEN
856  number_of_vertices=3
857  ELSE
858  number_of_vertices=4
859  ENDIF
860  ALLOCATE(gausspoints(number_of_vertices,max_gauss),stat=err)
861  IF(err/=0) CALL flagerror("Could not allocate weights",err,error,*999)
862  ALLOCATE(gaussweights(max_gauss),stat=err)
863  IF(err/=0) CALL flagerror("Could not allocate weights",err,error,*999)
864 
865  CALL gauss_simplex(order,number_of_vertices,numbergausspoints,gausspoints,gaussweights,err,error,*999)
866  CASE DEFAULT
867  local_error="Basis type "// &
868  & trim(number_to_vstring(basis%TYPE,"*",err,error))// &
869  & " is invalid or not implemented"
870  CALL flagerror(local_error,err,error,*999)
871  END SELECT
872  ELSE
873  CALL flagerror("Basis is not associated",err,error,*999)
874  ENDIF
875 
876  exits("BASIS_GAUSS_POINTS_CALCULATE_DP")
877  RETURN
878 999 errorsexits("BASIS_GAUSS_POINTS_CALCULATE_DP",err,error)
879  RETURN
880 
881  END SUBROUTINE basis_gauss_points_calculate_dp
882 
883  !
884  !================================================================================================================================
885  !
886 
889  RECURSIVE SUBROUTINE basis_family_destroy(USER_NUMBER,FAMILY_NUMBER,ERR,ERROR,*)
890 
891  !Argument variables
892  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
893  INTEGER(INTG), INTENT(IN) :: FAMILY_NUMBER
894  INTEGER(INTG), INTENT(OUT) :: ERR
895  TYPE(varying_string), INTENT(OUT) :: ERROR
896  !Local Variables
897  INTEGER(INTG) :: count,nb
898  TYPE(basis_type), POINTER :: BASIS
899  TYPE(basis_ptr_type), POINTER :: NEW_SUB_BASES(:)
900 
901  enters("BASIS_FAMILY_DESTROY",err,error,*999)
902 
903  NULLIFY(basis)
904  CALL basis_family_number_find(user_number,family_number,basis,err,error,*999)
905  IF(ASSOCIATED(basis)) THEN
906 
907 !!NOTE: We have to find a pointer to the basis to destroy within this routine rather than passing in a pointer to a
908 !!DESTROY_BASIS_PTR type routine because we need to change BASIS%SUB_BASES of the PARENT basis and this would violate section
909 !!12.4.1.6 of the Fortran standard if the dummy BASIS pointer argument was associated with the SUB_BASES(x)%PTR actual
910 !!argument.
911 
912  IF(basis%NUMBER_OF_SUB_BASES==0) THEN
913  !No more sub-bases so delete this instance
914  IF(ASSOCIATED(basis%PARENT_BASIS)) THEN
915  !Sub-basis function - delete this instance from the PARENT_BASIS
916  NULLIFY(new_sub_bases)
917  IF(basis%PARENT_BASIS%NUMBER_OF_SUB_BASES>1) THEN
918  !If the parent basis has more than one sub basis then remove this instance from its sub-bases list
919  ALLOCATE(new_sub_bases(basis%PARENT_BASIS%NUMBER_OF_SUB_BASES-1),stat=err)
920  IF(err/=0) CALL flagerror("Could not allocate new sub-bases",err,error,*999)
921  count=0
922  DO nb=1,basis%PARENT_BASIS%NUMBER_OF_SUB_BASES
923  IF(basis%PARENT_BASIS%SUB_BASES(nb)%PTR%USER_NUMBER==basis%USER_NUMBER.AND. &
924  & basis%PARENT_BASIS%SUB_BASES(nb)%PTR%FAMILY_NUMBER/=basis%FAMILY_NUMBER) THEN
925  count=count+1
926  new_sub_bases(count)%PTR=>basis%PARENT_BASIS%SUB_BASES(nb)%PTR
927  ENDIF
928  ENDDO
929  ENDIF
930  basis%PARENT_BASIS%NUMBER_OF_SUB_BASES=basis%PARENT_BASIS%NUMBER_OF_SUB_BASES-1
931  IF(ASSOCIATED(basis%PARENT_BASIS%SUB_BASES)) DEALLOCATE(basis%PARENT_BASIS%SUB_BASES)
932  basis%PARENT_BASIS%SUB_BASES=>new_sub_bases
933  ELSE
934  !Master basis function - delete this instance from BASIS_FUNCTIONS
935  NULLIFY(new_sub_bases)
936  IF(basis_functions%NUMBER_BASIS_FUNCTIONS>1) THEN
937  !If there is more than one basis defined then remove this instance from the basis functions
938  ALLOCATE(new_sub_bases(basis_functions%NUMBER_BASIS_FUNCTIONS-1),stat=err)
939  IF(err/=0) CALL flagerror("Could not allocate new sub-bases",err,error,*999)
940  count=0
941  DO nb=1,basis_functions%NUMBER_BASIS_FUNCTIONS
942  IF(basis_functions%BASES(nb)%PTR%USER_NUMBER/=basis%USER_NUMBER.AND. &
943  & basis_functions%BASES(nb)%PTR%FAMILY_NUMBER==0) THEN
944  count=count+1
945  new_sub_bases(count)%PTR=>basis_functions%BASES(nb)%PTR
946  ENDIF
947  ENDDO
948  ENDIF
949  IF(ASSOCIATED(basis_functions%BASES)) DEALLOCATE(basis_functions%BASES)
950  basis_functions%NUMBER_BASIS_FUNCTIONS=basis_functions%NUMBER_BASIS_FUNCTIONS-1
951  basis_functions%BASES=>new_sub_bases
952  ENDIF
953 
954  CALL basis_finalise(basis,err,error,*999)
955 
956  ELSE
957  !Recursively delete sub-bases first
958  DO WHILE(basis%NUMBER_OF_SUB_BASES>0)
959  CALL basis_family_destroy(basis%SUB_BASES(1)%PTR%USER_NUMBER,basis%SUB_BASES(1)%PTR%FAMILY_NUMBER,err,error,*999)
960  ENDDO
961  !Now delete this instance
962  CALL basis_family_destroy(user_number,family_number,err,error,*999)
963  ENDIF
964 
965  ELSE
966  CALL flagerror("Basis user number does not exist",err,error,*999)
967  ENDIF
968 
969  exits("BASIS_FAMILY_DESTROY")
970  RETURN
971 999 errorsexits("BASIS_FAMILY_DESTROY",err,error)
972  RETURN 1
973  END SUBROUTINE basis_family_destroy
974 
975  !
976  !================================================================================================================================
977  !
978 
981  RECURSIVE SUBROUTINE basis_family_number_find(USER_NUMBER,FAMILY_NUMBER,BASIS,ERR,ERROR,*)
982 
983  !Argument variables
984  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
985  INTEGER(INTG), INTENT(IN) :: FAMILY_NUMBER
986  TYPE(basis_type), POINTER :: BASIS
987  INTEGER(INTG), INTENT(OUT) :: ERR
988  TYPE(varying_string), INTENT(OUT) :: ERROR
989  !Local Variables
990  INTEGER(INTG) :: nb,nsb
991  TYPE(basis_type), POINTER :: SUB_BASIS
992 
993  enters("BASIS_FAMILY_NUMBER_FIND",err,error,*999)
994 
995  NULLIFY(basis)
996  nb=1
997  DO WHILE(nb<=basis_functions%NUMBER_BASIS_FUNCTIONS.AND..NOT.ASSOCIATED(basis))
998  IF(basis_functions%BASES(nb)%PTR%USER_NUMBER==user_number) THEN
999  IF(family_number==0) THEN
1000  basis=>basis_functions%BASES(nb)%PTR
1001  ELSE
1002 !!TODO: \todo This only works for one level of sub-bases at the moment
1003  nsb=1
1004  DO WHILE(nsb<=basis_functions%BASES(nb)%PTR%NUMBER_OF_SUB_BASES.AND..NOT.ASSOCIATED(basis))
1005  sub_basis=>basis_functions%BASES(nb)%PTR%SUB_BASES(nsb)%PTR
1006  IF(sub_basis%FAMILY_NUMBER==family_number) THEN
1007  basis=>sub_basis
1008  ELSE
1009  nsb=nsb+1
1010  ENDIF
1011  ENDDO
1012  ENDIF
1013  ELSE
1014  nb=nb+1
1015  ENDIF
1016  END DO
1017 
1018  exits("BASIS_FAMILY_NUMBER_FIND")
1019  RETURN
1020 999 errorsexits("BASIS_FAMILY_NUMBER_FIND",err,error)
1021  RETURN 1
1022  END SUBROUTINE basis_family_number_find
1023 
1024  !
1025  !================================================================================================================================
1026  !
1027 
1029  SUBROUTINE basis_finalise(BASIS,ERR,ERROR,*)
1030 
1031  !Argument variables
1032  TYPE(basis_type), POINTER :: BASIS
1033  INTEGER(INTG), INTENT(OUT) :: ERR
1034  TYPE(varying_string), INTENT(OUT) :: ERROR
1035  !Local Variables
1036 
1037  enters("BASIS_FINALISE",err,error,*999)
1038 
1039  IF(ASSOCIATED(basis)) THEN
1040  IF(ALLOCATED(basis%INTERPOLATION_XI)) DEALLOCATE(basis%INTERPOLATION_XI)
1041  IF(ALLOCATED(basis%INTERPOLATION_TYPE)) DEALLOCATE(basis%INTERPOLATION_TYPE)
1042  IF(ALLOCATED(basis%INTERPOLATION_ORDER)) DEALLOCATE(basis%INTERPOLATION_ORDER)
1043  IF(ALLOCATED(basis%COLLAPSED_XI)) DEALLOCATE(basis%COLLAPSED_XI)
1044  IF(ALLOCATED(basis%NODE_AT_COLLAPSE)) DEALLOCATE(basis%NODE_AT_COLLAPSE)
1045  CALL basis_quadrature_finalise(basis,err,error,*999)
1046  IF(ALLOCATED(basis%NUMBER_OF_NODES_XIC)) DEALLOCATE(basis%NUMBER_OF_NODES_XIC)
1047  IF(ALLOCATED(basis%NUMBER_OF_DERIVATIVES)) DEALLOCATE(basis%NUMBER_OF_DERIVATIVES)
1048  IF(ALLOCATED(basis%NODE_POSITION_INDEX)) DEALLOCATE(basis%NODE_POSITION_INDEX)
1049  IF(ALLOCATED(basis%NODE_POSITION_INDEX_INV)) DEALLOCATE(basis%NODE_POSITION_INDEX_INV)
1050  IF(ALLOCATED(basis%DERIVATIVE_ORDER_INDEX)) DEALLOCATE(basis%DERIVATIVE_ORDER_INDEX)
1051  IF(ALLOCATED(basis%DERIVATIVE_ORDER_INDEX_INV)) DEALLOCATE(basis%DERIVATIVE_ORDER_INDEX_INV)
1052  IF(ALLOCATED(basis%PARTIAL_DERIVATIVE_INDEX)) DEALLOCATE(basis%PARTIAL_DERIVATIVE_INDEX)
1053  IF(ALLOCATED(basis%ELEMENT_PARAMETER_INDEX)) DEALLOCATE(basis%ELEMENT_PARAMETER_INDEX)
1054  IF(ALLOCATED(basis%ELEMENT_PARAMETER_INDEX_INV)) DEALLOCATE(basis%ELEMENT_PARAMETER_INDEX_INV)
1055  IF(ALLOCATED(basis%LOCAL_LINE_XI_DIRECTION)) DEALLOCATE(basis%LOCAL_LINE_XI_DIRECTION)
1056  IF(ALLOCATED(basis%NUMBER_OF_NODES_IN_LOCAL_LINE)) DEALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE)
1057  IF(ALLOCATED(basis%NODE_NUMBERS_IN_LOCAL_LINE)) DEALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE)
1058  IF(ALLOCATED(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE)) DEALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE)
1059  IF(ALLOCATED(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE)) DEALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE)
1060  IF(ALLOCATED(basis%LOCAL_FACE_XI_DIRECTION)) DEALLOCATE(basis%LOCAL_FACE_XI_DIRECTION)
1061  IF(ALLOCATED(basis%NUMBER_OF_NODES_IN_LOCAL_FACE)) DEALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_FACE)
1062  IF(ALLOCATED(basis%NODE_NUMBERS_IN_LOCAL_FACE)) DEALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_FACE)
1063  IF(ALLOCATED(basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE)) DEALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE)
1064  IF(ALLOCATED(basis%ELEMENT_PARAMETERS_IN_LOCAL_FACE)) DEALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_FACE)
1065  IF(ALLOCATED(basis%LOCAL_XI_NORMAL)) DEALLOCATE(basis%LOCAL_XI_NORMAL)
1066  IF(ASSOCIATED(basis%LINE_BASES)) DEALLOCATE(basis%LINE_BASES)
1067  IF(ASSOCIATED(basis%FACE_BASES)) DEALLOCATE(basis%FACE_BASES)
1068  IF(ASSOCIATED(basis%SUB_BASES)) DEALLOCATE(basis%SUB_BASES)
1069  DEALLOCATE(basis)
1070  ENDIF
1071 
1072  exits("BASIS_FINALISE")
1073  RETURN
1074 999 errorsexits("BASIS_FINALISE",err,error)
1075  RETURN 1
1076 
1077  END SUBROUTINE basis_finalise
1078 
1079  !
1080  !================================================================================================================================
1081  !
1082 
1084  SUBROUTINE basis_initialise(BASIS,ERR,ERROR,*)
1085 
1086  !Argument variables
1087  TYPE(basis_type), POINTER :: BASIS
1088  INTEGER(INTG), INTENT(OUT) :: ERR
1089  TYPE(varying_string), INTENT(OUT) :: ERROR
1090  !Local Variables
1091 
1092  enters("BASIS_INITIALISE",err,error,*999)
1093 
1094  IF(ASSOCIATED(basis)) THEN
1095  basis%USER_NUMBER=0
1096  basis%GLOBAL_NUMBER=0
1097  basis%FAMILY_NUMBER=0
1098  basis%BASIS_FINISHED=.false.
1099  basis%HERMITE=.false.
1100  basis%TYPE=0
1101  basis%NUMBER_OF_XI=0
1102  basis%NUMBER_OF_XI_COORDINATES=0
1103  basis%DEGENERATE=.false.
1104  basis%NUMBER_OF_COLLAPSED_XI=0
1105  basis%NUMBER_OF_PARTIAL_DERIVATIVES=0
1106  basis%NUMBER_OF_NODES=0
1107  basis%NUMBER_OF_ELEMENT_PARAMETERS=0
1108  basis%MAXIMUM_NUMBER_OF_DERIVATIVES=0
1109  basis%NUMBER_OF_LOCAL_LINES=0
1110  basis%NUMBER_OF_LOCAL_FACES=0
1111  NULLIFY(basis%LINE_BASES)
1112  NULLIFY(basis%FACE_BASES)
1113  basis%NUMBER_OF_SUB_BASES=0
1114  NULLIFY(basis%SUB_BASES)
1115  NULLIFY(basis%PARENT_BASIS)
1116  ELSE
1117  CALL flagerror("Basis is not associated.",err,error,*999)
1118  ENDIF
1119 
1120  exits("BASIS_INITIALISE")
1121  RETURN
1122 999 errorsexits("BASIS_INITIALISE",err,error)
1123  RETURN 1
1124 
1125  END SUBROUTINE basis_initialise
1126 
1127  !
1128  !================================================================================================================================
1129  !
1130 
1134  FUNCTION basis_interpolate_gauss_dp(BASIS,PARTIAL_DERIV_INDEX,QUADRATURE_SCHEME,GAUSS_POINT_NUMBER,ELEMENT_PARAMETERS,ERR,ERROR)
1135 
1136  !Argument variables
1137  TYPE(basis_type), POINTER :: BASIS
1138  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIV_INDEX
1139  INTEGER(INTG), INTENT(IN) :: QUADRATURE_SCHEME
1140  INTEGER(INTG), INTENT(IN) :: GAUSS_POINT_NUMBER
1141  REAL(DP), INTENT(IN) :: ELEMENT_PARAMETERS(:)
1142  INTEGER(INTG), INTENT(OUT) :: ERR
1143  TYPE(varying_string), INTENT(OUT) :: ERROR
1144  !Function variable
1145  REAL(DP) :: BASIS_INTERPOLATE_GAUSS_DP
1146  !Local Variables
1147  INTEGER(INTG) :: ns
1148  TYPE(quadrature_scheme_type), POINTER :: BASIS_QUADRATURE_SCHEME
1149  TYPE(varying_string) :: LOCAL_ERROR
1150 
1151  enters("BASIS_INTERPOLATE_GAUSS_DP",err,error,*999)
1152 
1153  basis_interpolate_gauss_dp=0.0_dp
1154  IF(ASSOCIATED(basis)) THEN
1155  IF(quadrature_scheme>0.AND.quadrature_scheme<=basis_number_of_quadrature_scheme_types) THEN
1156  basis_quadrature_scheme=>basis%QUADRATURE%QUADRATURE_SCHEME_MAP(quadrature_scheme)%PTR
1157  IF(ASSOCIATED(basis_quadrature_scheme)) THEN
1158  IF(gauss_point_number>0.AND.gauss_point_number<=basis_quadrature_scheme%NUMBER_OF_GAUSS) THEN
1159  IF(partial_deriv_index>0.AND.partial_deriv_index<=basis%NUMBER_OF_PARTIAL_DERIVATIVES) THEN
1160  DO ns=1,basis%NUMBER_OF_ELEMENT_PARAMETERS
1161  basis_interpolate_gauss_dp=basis_interpolate_gauss_dp+ &
1162  & basis_quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_deriv_index,gauss_point_number)* &
1163  & element_parameters(ns)
1164  ENDDO !ns
1165  ELSE
1166  local_error="The partial derivative index of "//trim(number_to_vstring(partial_deriv_index,"*",err,error))// &
1167  & " is invalid. It must be between 1 and "// &
1168  & trim(number_to_vstring(basis%NUMBER_OF_PARTIAL_DERIVATIVES,"*",err,error))
1169  CALL flagerror(local_error,err,error,*999)
1170  ENDIF
1171  ENDIF
1172  ELSE
1173  CALL flagerror("The quadrature scheme has not been created",err,error,*999)
1174  ENDIF
1175  ELSE
1176  local_error="The quadrature scheme type number of "//trim(number_to_vstring(quadrature_scheme,"*",err,error))// &
1177  & " is invalid. It must be between 1 and "// &
1178  & trim(number_to_vstring(basis_number_of_quadrature_scheme_types,"*",err,error))
1179  CALL flagerror(local_error,err,error,*999)
1180  ENDIF
1181  ELSE
1182  CALL flagerror("Basis is not associated",err,error,*999)
1183  ENDIF
1184 
1185  exits("BASIS_INTERPOLATE_GAUSS_DP")
1186  RETURN
1187 999 errorsexits("BASIS_INTERPOLATE_GAUSS_DP",err,error)
1188  END FUNCTION basis_interpolate_gauss_dp
1189 
1190  !
1191  !================================================================================================================================
1192  !
1193 
1197  FUNCTION basis_interpolate_local_face_gauss_dp(BASIS,PARTIAL_DERIV_INDEX,QUADRATURE_SCHEME, &
1198  & local_face_number,gauss_point_number,face_parameters,err,error)
1199 
1200  !Argument variables
1201  TYPE(basis_type), POINTER :: BASIS
1202  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIV_INDEX
1203  INTEGER(INTG), INTENT(IN) :: QUADRATURE_SCHEME
1204  INTEGER(INTG), INTENT(IN) :: LOCAL_FACE_NUMBER
1205  INTEGER(INTG), INTENT(IN) :: GAUSS_POINT_NUMBER
1206  REAL(DP), INTENT(IN) :: FACE_PARAMETERS(:)
1207  INTEGER(INTG), INTENT(OUT) :: ERR
1208  TYPE(varying_string), INTENT(OUT) :: ERROR
1209  !Function variable
1210  REAL(DP) :: BASIS_INTERPOLATE_LOCAL_FACE_GAUSS_DP
1211  !Local Variables
1212  INTEGER(INTG) :: ns
1213  TYPE(quadrature_scheme_type), POINTER :: BASIS_QUADRATURE_SCHEME
1214  TYPE(varying_string) :: LOCAL_ERROR
1215 
1216  enters("BASIS_INTERPOLATE_LOCAL_FACE_GAUSS_DP",err,error,*999)
1217 
1218  basis_interpolate_local_face_gauss_dp=0.0_dp
1219  IF(ASSOCIATED(basis)) THEN
1220  IF(quadrature_scheme>0.AND.quadrature_scheme<=basis_number_of_quadrature_scheme_types) THEN
1221  basis_quadrature_scheme=>basis%QUADRATURE%QUADRATURE_SCHEME_MAP(quadrature_scheme)%PTR
1222  IF(ASSOCIATED(basis_quadrature_scheme)) THEN
1223  IF(basis%QUADRATURE%EVALUATE_FACE_GAUSS) THEN !Alternartively, can check whether scheme's face arrays are allocated?
1224  IF(local_face_number>0.AND.local_face_number<=basis%NUMBER_OF_LOCAL_FACES) THEN
1225  IF(gauss_point_number>0.AND.gauss_point_number<=basis_quadrature_scheme%NUMBER_OF_FACE_GAUSS(local_face_number)) THEN
1226  IF(partial_deriv_index>0.AND.partial_deriv_index<=basis%NUMBER_OF_PARTIAL_DERIVATIVES) THEN
1227  DO ns=1,basis%NUMBER_OF_ELEMENT_PARAMETERS
1228  basis_interpolate_local_face_gauss_dp=basis_interpolate_local_face_gauss_dp+ &
1229  & basis_quadrature_scheme%FACE_GAUSS_BASIS_FNS(ns,partial_deriv_index,gauss_point_number,local_face_number)* &
1230  & face_parameters(ns)
1231  ENDDO !ns
1232  ELSE
1233  local_error="The partial derivative index of "//trim(number_to_vstring(partial_deriv_index,"*",err,error))// &
1234  & " is invalid. It must be between 1 and "// &
1235  & trim(number_to_vstring(basis%NUMBER_OF_PARTIAL_DERIVATIVES,"*",err,error))
1236  CALL flagerror(local_error,err,error,*999)
1237  ENDIF
1238  ENDIF
1239  ELSE
1240  CALL flagerror("The local face number index is invalid.",err,error,*999)
1241  ENDIF
1242  ELSE
1243  CALL flagerror("The face gauss interpolation scheme has not been created",err,error,*999)
1244  ENDIF
1245  ELSE
1246  CALL flagerror("The quadrature scheme has not been created",err,error,*999)
1247  ENDIF
1248  ELSE
1249  local_error="The quadrature scheme type number of "//trim(number_to_vstring(quadrature_scheme,"*",err,error))// &
1250  & " is invalid. It must be between 1 and "// &
1251  & trim(number_to_vstring(basis_number_of_quadrature_scheme_types,"*",err,error))
1252  CALL flagerror(local_error,err,error,*999)
1253  ENDIF
1254  ELSE
1255  CALL flagerror("Basis is not associated",err,error,*999)
1256  ENDIF
1257 
1258  exits("BASIS_INTERPOLATE_LOCAL_FACE_GAUSS_DP")
1259  RETURN
1260 999 errorsexits("BASIS_INTERPOLATE_LOCAL_FACE_GAUSS_DP",err,error)
1261  END FUNCTION basis_interpolate_local_face_gauss_dp
1262 
1263  !
1264  !================================================================================================================================
1265  !
1266 
1271  FUNCTION basis_interpolate_xi_dp(BASIS,PARTIAL_DERIV_INDEX,XI,ELEMENT_PARAMETERS,ERR,ERROR)
1272 
1273  !Argument variables
1274  TYPE(basis_type), POINTER :: BASIS
1275  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIV_INDEX
1276  REAL(DP), INTENT(IN) :: XI(:)
1277  REAL(DP), INTENT(IN) :: ELEMENT_PARAMETERS(:)
1278  INTEGER(INTG), INTENT(OUT) :: ERR
1279  TYPE(varying_string), INTENT(OUT) :: ERROR
1280  !Function variable
1281  REAL(DP) :: BASIS_INTERPOLATE_XI_DP
1282  !Local Variables
1283  INTEGER(INTG) :: nn,nk,ns
1284  REAL(DP) :: XIL(size(xi,1)+1)
1285  TYPE(varying_string) :: LOCAL_ERROR
1286 
1287  enters("BASIS_INTERPOLATE_XI_DP",err,error,*999)
1288 
1289  basis_interpolate_xi_dp=0.0_dp
1290  IF(ASSOCIATED(basis)) THEN
1291  SELECT CASE(basis%TYPE)
1292  CASE(basis_lagrange_hermite_tp_type)
1293  ns=0
1294  DO nn=1,basis%NUMBER_OF_NODES
1295  DO nk=1,basis%NUMBER_OF_DERIVATIVES(nn)
1296  ns=ns+1
1297  basis_interpolate_xi_dp=basis_interpolate_xi_dp+ &
1298  & basis_lhtp_basis_evaluate(basis,nn,nk,partial_deriv_index,xi,err,error)* &
1299  & element_parameters(ns)
1300  ENDDO !nk
1301  ENDDO !nn
1302  IF(err/=0) GOTO 999
1303  CASE(basis_simplex_type)
1304  !Create the area coordinates from the xi coordinates
1305  xil(1:SIZE(xi,1))=1.0_dp-xi
1306  xil(SIZE(xi,1)+1)=sum(xi)-(SIZE(xi,1)-1.0_dp)
1307  ns=0
1308  DO nn=1,basis%NUMBER_OF_NODES
1309  ns=ns+1
1310  basis_interpolate_xi_dp=basis_interpolate_xi_dp+ &
1311  & basis_simplex_basis_evaluate(basis,nn,partial_deriv_index,xil,err,error)* &
1312  & element_parameters(ns)
1313  ENDDO !nn
1314  IF(err/=0) GOTO 999
1315  CASE(basis_radial_type)
1316 
1317  IF(err/=0) GOTO 999
1318  CASE DEFAULT
1319  local_error="Basis type "//trim(number_to_vstring(basis%TYPE,"*",err,error))//" is invalid or not implemented"
1320  CALL flagerror(local_error,err,error,*999)
1321  END SELECT
1322  ELSE
1323  CALL flagerror("Basis is not associated",err,error,*999)
1324  ENDIF
1325 
1326  exits("BASIS_INTERPOLATE_XI_DP")
1327  RETURN
1328 999 errorsexits("BASIS_INTERPOLATE_XI_DP",err,error)
1329  END FUNCTION basis_interpolate_xi_dp
1330 
1331  !
1332  !================================================================================================================================
1333  !
1334 
1336  SUBROUTINE basis_interpolation_xi_get(BASIS,INTERPOLATION_XI,ERR,ERROR,*)
1337 
1338  !Argument variables
1339  TYPE(basis_type), POINTER :: BASIS
1340  INTEGER(INTG), INTENT(OUT) :: INTERPOLATION_XI(:)
1341  INTEGER(INTG), INTENT(OUT) :: ERR
1342  TYPE(varying_string), INTENT(OUT) :: ERROR
1343  !Local Variables
1344  TYPE(varying_string) :: LOCAL_ERROR
1345 
1346  enters("BASIS_INTERPOLATION_XI_GET",err,error,*999)
1347 
1348  IF(ASSOCIATED(basis)) THEN
1349  IF(basis%BASIS_FINISHED) THEN
1350  IF(SIZE(interpolation_xi,1)>=SIZE(basis%INTERPOLATION_XI,1)) THEN
1351  interpolation_xi=basis%INTERPOLATION_XI
1352  ELSE
1353  local_error="The size of INTERPOLATION_XI is too small. The supplied size is "// &
1354  & trim(number_to_vstring(SIZE(interpolation_xi,1),"*",err,error))//" and it needs to be >= "// &
1355  & trim(number_to_vstring(SIZE(basis%INTERPOLATION_XI,1),"*",err,error))//"."
1356  CALL flagerror(local_error,err,error,*999)
1357  ENDIF
1358  ELSE
1359  CALL flagerror("Basis has not been finished.",err,error,*999)
1360  ENDIF
1361  ELSE
1362  CALL flagerror("Basis is not associated.",err,error,*999)
1363  ENDIF
1364 
1365  exits("BASIS_INTERPOLATION_XI_GET")
1366  RETURN
1367 999 errorsexits("BASIS_INTERPOLATION_XI_GET",err,error)
1368  RETURN
1369  END SUBROUTINE basis_interpolation_xi_get
1370 
1371 
1372  !
1373  !================================================================================================================================
1374  !
1375 
1377  SUBROUTINE basis_interpolation_xi_set_number(USER_NUMBER,INTERPOLATION_XI,ERR,ERROR,*)
1378 
1379  !Argument variables
1380  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
1381  INTEGER(INTG), INTENT(IN) :: INTERPOLATION_XI(:)
1382  INTEGER(INTG), INTENT(OUT) :: ERR
1383  TYPE(varying_string), INTENT(OUT) :: ERROR
1384  !Local Variables
1385  TYPE(basis_type), POINTER :: BASIS
1386 
1387  enters("BASIS_INTERPOLATION_XI_SET_NUMBER",err,error,*999)
1388 
1389  CALL basis_user_number_find(user_number,basis,err,error,*999)
1390  CALL basis_interpolation_xi_set(basis,interpolation_xi,err,error,*999)
1391 
1392  exits("BASIS_INTERPOLATION_XI_SET_NUMBER")
1393  RETURN
1394 999 errorsexits("BASIS_INTERPOLATION_XI_SET_NUMBER",err,error)
1395  RETURN 1
1396  END SUBROUTINE basis_interpolation_xi_set_number
1397 
1398  !
1399  !================================================================================================================================
1400  !
1401 
1403  SUBROUTINE basis_interpolation_xi_set_ptr(BASIS,INTERPOLATION_XI,ERR,ERROR,*)
1404 
1405  !Argument variables
1406  TYPE(basis_type), POINTER :: BASIS
1407  INTEGER(INTG), INTENT(IN) :: INTERPOLATION_XI(:)
1408  INTEGER(INTG), INTENT(OUT) :: ERR
1409  TYPE(varying_string), INTENT(OUT) :: ERROR
1410  !Local Variables
1411  INTEGER(INTG) :: ni,LAST_INTERP
1412  TYPE(varying_string) :: LOCAL_ERROR
1413 
1414  enters("BASIS_INTERPOLATION_XI_SET_PTR",err,error,*999)
1415 
1416  IF(ASSOCIATED(basis)) THEN
1417  IF(basis%BASIS_FINISHED) THEN
1418  CALL flagerror("Basis has been finished",err,error,*999)
1419  ELSE
1420  IF(SIZE(interpolation_xi,1)==basis%NUMBER_OF_XI) THEN
1421  !Check the input values
1422  SELECT CASE(basis%TYPE)
1423  CASE(basis_lagrange_hermite_tp_type)
1424  DO ni=1,basis%NUMBER_OF_XI
1425  SELECT CASE(interpolation_xi(ni))
1426  CASE(basis_linear_lagrange_interpolation,basis_quadratic_lagrange_interpolation,basis_cubic_lagrange_interpolation, &
1427  & basis_cubic_hermite_interpolation,basis_quadratic1_hermite_interpolation,basis_quadratic2_hermite_interpolation)
1428  !Do nothing
1429  CASE DEFAULT
1430  local_error="Interpolation xi value "//trim(number_to_vstring(interpolation_xi(ni),"*",err,error))// &
1431  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid for a Lagrange-Hermite TP basis."
1432  CALL flagerror(local_error,err,error,*999)
1433  END SELECT
1434  ENDDO !ni
1435  CASE(basis_simplex_type)
1436  last_interp=interpolation_xi(1)
1437  DO ni=1,basis%NUMBER_OF_XI
1438  SELECT CASE(interpolation_xi(ni))
1439  CASE(basis_linear_simplex_interpolation,basis_quadratic_simplex_interpolation,basis_cubic_simplex_interpolation)
1440  IF(interpolation_xi(ni)/=last_interp) THEN
1441  CALL flagerror("The interpolation xi value must be the same for all xi directions for a simplex basis.", &
1442  & err,error,*999)
1443  ENDIF
1444  CASE DEFAULT
1445  local_error="Interpolation xi value "//trim(number_to_vstring(interpolation_xi(ni),"*",err,error))// &
1446  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid for a simplex basis."
1447  CALL flagerror(local_error,err,error,*999)
1448  END SELECT
1449  ENDDO !ni
1450  CASE DEFAULT
1451  CALL flagerror("Invalid basis type or not implemented",err,error,*999)
1452  END SELECT
1453  !Set the interpolation xi
1454  basis%INTERPOLATION_XI(1:basis%NUMBER_OF_XI)=interpolation_xi(1:basis%NUMBER_OF_XI)
1455  ELSE
1456  local_error="The size of the interpolation xi array ("// &
1457  & trim(number_to_vstring(SIZE(interpolation_xi,1),"*",err,error))//") does not match the number of xi directions ("// &
1458  & trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))//") for basis number "// &
1459  & trim(number_to_vstring(basis%USER_NUMBER,"*",err,error))//"."
1460  CALL flagerror(local_error,err,error,*999)
1461  ENDIF
1462  ENDIF
1463  ELSE
1464  CALL flagerror("Basis is not associated.",err,error,*999)
1465  ENDIF
1466 
1467  exits("BASIS_INTERPOLATION_XI_SET_PTR")
1468  RETURN
1469 999 errorsexits("BASIS_INTERPOLATION_XI_SET_PTR",err,error)
1470  RETURN 1
1471  END SUBROUTINE basis_interpolation_xi_set_ptr
1472 
1473  !
1474  !================================================================================================================================
1475  !
1476 
1479  SUBROUTINE basis_lhtpbasiscreate(basis,err,error,*)
1480 
1481  !Argument variables
1482  TYPE(basis_type), POINTER :: basis
1483  INTEGER(INTG), INTENT(OUT) :: err
1484  TYPE(varying_string), INTENT(OUT) :: error
1485  !Local Variables
1486  INTEGER(INTG) :: maximumNumberOfNodes,numberOfDerivatives,xiIdx,xiIdx1,xiIdx2,xiIdx3,derivativeIdx,localNode,localLineNodeIdx, &
1487  & elementParameter,oldNumberOfDerivatives,position(4),collapsedPosition(3),maximumNodeExtent(3),collapsedXi(3), &
1488  & numberOfNodes,numberOfLocalLines,nodeCount,specialNodeCount,nodesInLine(4),numberOfLocalFaces,localFaceIdx, &
1489  & localNodeIdx,localNodeIdx1,localNodeIdx2,localNodeIdx3,directionIdx,localFaceDerivative,localNodeCount, &
1490  & localLineParameter,localFaceParameter
1491  LOGICAL, ALLOCATABLE :: nodeAtCollapse(:)
1492  LOGICAL :: atCollapse,collapsedFace,firstCollapsedPosition
1493 
1494  enters("Basis_LHTPBasisCreate",err,error,*999)
1495 
1496  IF(ASSOCIATED(basis)) THEN
1497  IF(basis%type==basis_lagrange_hermite_tp_type) THEN
1498  basis%NUMBER_OF_XI_COORDINATES=basis%NUMBER_OF_XI
1499  basis%NUMBER_OF_PARTIAL_DERIVATIVES=basis%NUMBER_OF_XI_COORDINATES**2+2
1500  ALLOCATE(basis%INTERPOLATION_TYPE(basis%NUMBER_OF_XI_COORDINATES),stat=err)
1501  IF(err/=0) CALL flagerror("Could not allocate basis interpolation type.",err,error,*999)
1502  ALLOCATE(basis%INTERPOLATION_ORDER(basis%NUMBER_OF_XI_COORDINATES),stat=err)
1503  IF(err/=0) CALL flagerror("Could not allocate basis interpolation order.",err,error,*999)
1504  ALLOCATE(basis%NUMBER_OF_NODES_XIC(basis%NUMBER_OF_XI_COORDINATES),stat=err)
1505  IF(err/=0) CALL flagerror("Could not allocate basis number of nodes xic.",err,error,*999)
1506  numberofnodes=1
1507  maximumnumberofnodes=0
1508  basis%degenerate=.false.
1509  basis%NUMBER_OF_COLLAPSED_XI=0
1510  DO xiidx=1,basis%NUMBER_OF_XI
1511  !Set up the interpolation types, orders and number of nodes in each xi from the user specified interpolation xi.
1512  SELECT CASE(basis%INTERPOLATION_XI(xiidx))
1513  CASE(basis_linear_lagrange_interpolation)
1514  basis%INTERPOLATION_TYPE(xiidx)=basis_lagrange_interpolation
1515  basis%INTERPOLATION_ORDER(xiidx)=basis_linear_interpolation_order
1516  basis%NUMBER_OF_NODES_XIC(xiidx)=2
1517  CASE(basis_quadratic_lagrange_interpolation)
1518  basis%INTERPOLATION_TYPE(xiidx)=basis_lagrange_interpolation
1519  basis%INTERPOLATION_ORDER(xiidx)=basis_quadratic_interpolation_order
1520  basis%NUMBER_OF_NODES_XIC(xiidx)=3
1521  CASE(basis_cubic_lagrange_interpolation)
1522  basis%INTERPOLATION_TYPE(xiidx)=basis_lagrange_interpolation
1523  basis%INTERPOLATION_ORDER(xiidx)=basis_cubic_interpolation_order
1524  basis%NUMBER_OF_NODES_XIC(xiidx)=4
1525  CASE(basis_cubic_hermite_interpolation)
1526  basis%INTERPOLATION_TYPE(xiidx)=basis_hermite_interpolation
1527  basis%INTERPOLATION_ORDER(xiidx)=basis_cubic_interpolation_order
1528  basis%NUMBER_OF_NODES_XIC(xiidx)=2
1529  CASE(basis_quadratic1_hermite_interpolation)
1530  basis%INTERPOLATION_TYPE(xiidx)=basis_hermite_interpolation
1531  basis%INTERPOLATION_ORDER(xiidx)=basis_quadratic1_interpolation_order
1532  basis%NUMBER_OF_NODES_XIC(xiidx)=2
1533  CASE(basis_quadratic2_hermite_interpolation)
1534  basis%INTERPOLATION_TYPE(xiidx)=basis_hermite_interpolation
1535  basis%INTERPOLATION_ORDER(xiidx)=basis_quadratic2_interpolation_order
1536  basis%NUMBER_OF_NODES_XIC(xiidx)=2
1537  CASE DEFAULT
1538  CALL flagerror("Invalid interpolation type",err,error,*999)
1539  END SELECT
1540  IF(basis%COLLAPSED_XI(xiidx)==basis_xi_collapsed) THEN
1541  basis%NUMBER_OF_COLLAPSED_XI=basis%NUMBER_OF_COLLAPSED_XI+1
1542  collapsedxi(basis%NUMBER_OF_COLLAPSED_XI)=xiidx
1543  basis%degenerate=.true.
1544  ENDIF
1545  numberofnodes=numberofnodes*basis%NUMBER_OF_NODES_XIC(xiidx)
1546  IF(basis%NUMBER_OF_NODES_XIC(xiidx)>maximumnumberofnodes) maximumnumberofnodes=basis%NUMBER_OF_NODES_XIC(xiidx)
1547  ENDDO !xiIdx
1548  !If a degenerate (collapsed) basis recalculate the number of nodes from the maximum posible number of nodes
1549  IF(basis%degenerate) THEN
1550  !Calculate the NODE_AT_COLLAPSE array.
1551  ALLOCATE(nodeatcollapse(numberofnodes),stat=err)
1552  IF(err/=0) CALL flagerror("Could not allocate at collapse",err,error,*999)
1553  position=1
1554  basis%NUMBER_OF_NODES=0
1555  !Loop over the maximum number of nodes which is currently set for the basis
1556  DO localnodeidx=1,numberofnodes
1557  atcollapse=.false.
1558  DO xiidx=1,basis%NUMBER_OF_XI
1559  IF(basis%COLLAPSED_XI(xiidx)==basis_collapsed_at_xi0.AND.position(xiidx)==1.OR. &
1560  & basis%COLLAPSED_XI(xiidx)==basis_collapsed_at_xi1.AND.position(xiidx)==basis%NUMBER_OF_NODES_XIC(xiidx)) &
1561  & THEN
1562  atcollapse=.true.
1563  firstcollapsedposition=all(position(collapsedxi(1:basis%NUMBER_OF_COLLAPSED_XI))==1)
1564  EXIT
1565  ENDIF
1566  ENDDO !xiIdx
1567  IF(atcollapse) THEN
1568  IF(firstcollapsedposition) THEN
1569  basis%NUMBER_OF_NODES=basis%NUMBER_OF_NODES+1
1570  nodeatcollapse(basis%NUMBER_OF_NODES)=.true.
1571  ENDIF
1572  ELSE
1573  basis%NUMBER_OF_NODES=basis%NUMBER_OF_NODES+1
1574  nodeatcollapse(basis%NUMBER_OF_NODES)=.false.
1575  ENDIF
1576  position(1)=position(1)+1
1577  DO xiidx=1,basis%NUMBER_OF_XI
1578  IF(position(xiidx)>basis%NUMBER_OF_NODES_XIC(xiidx)) THEN
1579  position(xiidx)=1
1580  position(xiidx+1)=position(xiidx+1)+1
1581  ENDIF
1582  ENDDO !xiIdx
1583  ENDDO !localNodeIdx
1584  CALL move_alloc(nodeatcollapse,basis%NODE_AT_COLLAPSE)
1585  ELSE
1586  basis%NUMBER_OF_NODES=numberofnodes
1587  ALLOCATE(basis%NODE_AT_COLLAPSE(basis%NUMBER_OF_NODES),stat=err)
1588  IF(err/=0) CALL flagerror("Could not allocate basis node at collapse.",err,error,*999)
1589  basis%NODE_AT_COLLAPSE=.false.
1590  collapsedxi(1)=1
1591  ENDIF
1592 
1593  ALLOCATE(basis%NODE_POSITION_INDEX(basis%NUMBER_OF_NODES,basis%NUMBER_OF_XI_COORDINATES),stat=err)
1594  IF(err/=0) CALL flagerror("Could not allocate basis node position index.",err,error,*999)
1595  SELECT CASE(basis%NUMBER_OF_XI_COORDINATES)
1596  CASE(1)
1597  ALLOCATE(basis%NODE_POSITION_INDEX_INV(maximumnumberofnodes,1,1,1),stat=err)
1598  CASE(2)
1599  ALLOCATE(basis%NODE_POSITION_INDEX_INV(maximumnumberofnodes,maximumnumberofnodes,1,1),stat=err)
1600  CASE(3)
1601  ALLOCATE(basis%NODE_POSITION_INDEX_INV(maximumnumberofnodes,maximumnumberofnodes,maximumnumberofnodes,1),stat=err)
1602  CASE DEFAULT
1603  CALL flagerror("Invalid number of xi coordinates.",err,error,*999)
1604  END SELECT
1605  IF(err/=0) CALL flagerror("Could not allocate node position index inverse.",err,error,*999)
1606  basis%NODE_POSITION_INDEX_INV=0
1607 
1608  !Determine the node position index and its inverse
1609  position=1
1610  collapsedposition=1
1611  localnode=0
1612  firstcollapsedposition=.true.
1613  DO localnodeidx1=1,numberofnodes
1614  atcollapse=.false.
1615  IF(basis%degenerate) THEN
1616  DO xiidx=1,basis%NUMBER_OF_XI
1617  IF(basis%COLLAPSED_XI(xiidx)==basis_collapsed_at_xi0.AND.position(xiidx)==1.OR. &
1618  & basis%COLLAPSED_XI(xiidx)==basis_collapsed_at_xi1.AND.position(xiidx)==basis%NUMBER_OF_NODES_XIC(xiidx)) &
1619  & THEN
1620  atcollapse=.true.
1621  firstcollapsedposition=all(position(collapsedxi(1:basis%NUMBER_OF_COLLAPSED_XI))==1)
1622  EXIT
1623  ENDIF
1624  ENDDO !xiIdx
1625  ENDIF
1626  IF(atcollapse) THEN
1627  IF(firstcollapsedposition) THEN
1628  localnode=localnode+1
1629  basis%NODE_POSITION_INDEX(localnode,1:basis%NUMBER_OF_XI)=position(1:basis%NUMBER_OF_XI)
1630  basis%NODE_POSITION_INDEX_INV(position(1),position(2),position(3),1)=localnode
1631  ELSE
1632  !The second node in the collapsed xi is set to the same node number as the first node in that xi direction.
1633  collapsedposition(1:basis%NUMBER_OF_XI)=position(1:basis%NUMBER_OF_XI)
1634  collapsedposition(collapsedxi(1:basis%NUMBER_OF_COLLAPSED_XI))=1
1635  basis%NODE_POSITION_INDEX_INV(position(1),position(2),position(3),1)= &
1636  & basis%NODE_POSITION_INDEX_INV(collapsedposition(1),collapsedposition(2),collapsedposition(3),1)
1637  ENDIF
1638  ELSE
1639  localnode=localnode+1
1640  basis%NODE_POSITION_INDEX(localnode,1:basis%NUMBER_OF_XI)=position(1:basis%NUMBER_OF_XI)
1641  basis%NODE_POSITION_INDEX_INV(position(1),position(2),position(3),1)=localnode
1642  ENDIF
1643  position(1)=position(1)+1
1644  DO xiidx=1,basis%NUMBER_OF_XI
1645  IF(position(xiidx)>basis%NUMBER_OF_NODES_XIC(xiidx)) THEN
1646  position(xiidx)=1
1647  position(xiidx+1)=position(xiidx+1)+1
1648  ENDIF
1649  ENDDO !xiIdx
1650  ENDDO !localNodeIdx1
1651  !Calculate the maximum number of derivatives and the number of element parameters
1652  basis%MAXIMUM_NUMBER_OF_DERIVATIVES=-1
1653  basis%NUMBER_OF_ELEMENT_PARAMETERS=0
1654  DO localnodeidx=1,basis%NUMBER_OF_NODES
1655  numberofderivatives=1
1656  DO xiidx=1,basis%NUMBER_OF_XI
1657  IF((.NOT.basis%NODE_AT_COLLAPSE(localnodeidx).OR.basis%COLLAPSED_XI(xiidx)==basis_not_collapsed).AND. &
1658  & basis%INTERPOLATION_TYPE(xiidx)==basis_hermite_interpolation.AND. &
1659  & (basis%INTERPOLATION_ORDER(xiidx)==basis_cubic_interpolation_order.OR. &
1660  & (basis%NODE_POSITION_INDEX(localnodeidx,xiidx)==1.AND. &
1661  & basis%INTERPOLATION_ORDER(xiidx)==basis_quadratic2_interpolation_order).OR. &
1662  & (basis%NODE_POSITION_INDEX(localnodeidx,xiidx)==2.AND. &
1663  & basis%INTERPOLATION_ORDER(xiidx)==basis_quadratic1_interpolation_order))) THEN
1664  !Derivative in this direction
1665  numberofderivatives=numberofderivatives*2
1666  ENDIF
1667  ENDDO !xiIdx
1668  basis%NUMBER_OF_ELEMENT_PARAMETERS=basis%NUMBER_OF_ELEMENT_PARAMETERS+numberofderivatives
1669  IF(numberofderivatives>basis%MAXIMUM_NUMBER_OF_DERIVATIVES) basis%MAXIMUM_NUMBER_OF_DERIVATIVES=numberofderivatives
1670  ENDDO !localNodeIdx
1671  !Now set up the number of derivatives and derivative order index
1672  ALLOCATE(basis%NUMBER_OF_DERIVATIVES(basis%NUMBER_OF_NODES),stat=err)
1673  IF(err/=0) CALL flagerror("Could not allocate number of derivatives.",err,error,*999)
1674  ALLOCATE(basis%DERIVATIVE_ORDER_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES, &
1675  & basis%NUMBER_OF_XI),stat=err)
1676  IF(err/=0) CALL flagerror("Could not allocate derivative order index.",err,error,*999)
1677  ALLOCATE(basis%DERIVATIVE_ORDER_INDEX_INV(first_part_deriv,first_part_deriv,first_part_deriv, &
1678  & basis%NUMBER_OF_NODES),stat=err)
1679  IF(err/=0) CALL flagerror("Could not allocate derivative order index inverse.",err,error,*999)
1680  ALLOCATE(basis%PARTIAL_DERIVATIVE_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES),stat=err)
1681  IF(err/=0) CALL flagerror("Could not allocate partial derivative index.",err,error,*999)
1682  ALLOCATE(basis%ELEMENT_PARAMETER_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES),stat=err)
1683  IF(err/=0) CALL flagerror("Could not allocate element parameter index.",err,error,*999)
1684  ALLOCATE(basis%ELEMENT_PARAMETER_INDEX_INV(2,basis%NUMBER_OF_ELEMENT_PARAMETERS),stat=err)
1685  IF(err/=0) CALL flagerror("Could not allocate element parameter index inverse.",err,error,*999)
1686  !Set the derivative order index and its inverse, the element parameter index and the partial derivative index.
1687  elementparameter=0
1688  basis%DERIVATIVE_ORDER_INDEX=0
1689  basis%DERIVATIVE_ORDER_INDEX_INV=0
1690  DO localnodeidx=1,basis%NUMBER_OF_NODES
1691  basis%NUMBER_OF_DERIVATIVES(localnodeidx)=1
1692  DO xiidx1=1,basis%NUMBER_OF_XI
1693  IF((.NOT.basis%NODE_AT_COLLAPSE(localnodeidx).OR.basis%COLLAPSED_XI(xiidx1)==basis_not_collapsed).AND. &
1694  & basis%INTERPOLATION_TYPE(xiidx1)==basis_hermite_interpolation.AND. &
1695  & (basis%INTERPOLATION_ORDER(xiidx1)==basis_cubic_interpolation_order.OR. &
1696  & (basis%NODE_POSITION_INDEX(localnodeidx,xiidx1)==1.AND. &
1697  & basis%INTERPOLATION_ORDER(xiidx1)==basis_quadratic2_interpolation_order).OR. &
1698  & (basis%NODE_POSITION_INDEX(localnodeidx,xiidx1)==2.AND. &
1699  & basis%INTERPOLATION_ORDER(xiidx1)==basis_quadratic1_interpolation_order))) THEN
1700  oldnumberofderivatives=basis%NUMBER_OF_DERIVATIVES(localnodeidx)
1701  basis%NUMBER_OF_DERIVATIVES(localnodeidx)=basis%NUMBER_OF_DERIVATIVES(localnodeidx)*2
1702  DO derivativeidx=1,oldnumberofderivatives
1703  basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,xiidx1)=no_part_deriv
1704  basis%DERIVATIVE_ORDER_INDEX(oldnumberofderivatives+derivativeidx,localnodeidx,xiidx1)=first_part_deriv
1705  DO xiidx2=1,xiidx1-1
1706  basis%DERIVATIVE_ORDER_INDEX(oldnumberofderivatives+derivativeidx,localnodeidx,xiidx2)= &
1707  & basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,xiidx2)
1708  ENDDO !xiIdx2
1709  ENDDO !derivativeIdx
1710  ELSE
1711  DO derivativeidx=1,basis%NUMBER_OF_DERIVATIVES(localnodeidx)
1712  basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,xiidx1)=no_part_deriv
1713  ENDDO !derivativeIdx
1714  ENDIF
1715  ENDDO !xiIdx1
1716 
1717  DO derivativeidx=1,basis%NUMBER_OF_DERIVATIVES(localnodeidx)
1718  elementparameter=elementparameter+1
1719  basis%ELEMENT_PARAMETER_INDEX(derivativeidx,localnodeidx)=elementparameter
1720  basis%ELEMENT_PARAMETER_INDEX_INV(1,elementparameter)=localnodeidx
1721  basis%ELEMENT_PARAMETER_INDEX_INV(2,elementparameter)=derivativeidx
1722  SELECT CASE(basis%NUMBER_OF_XI)
1723  CASE(1)
1724  basis%DERIVATIVE_ORDER_INDEX_INV(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,1),1,1,localnodeidx)= &
1725  & derivativeidx
1726  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,1))
1727  CASE(no_part_deriv)
1728  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=no_part_deriv
1729  CASE(first_part_deriv)
1730  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s1
1731  CASE DEFAULT
1732  CALL flagerror("Invalid derivative order index.",err,error,*999)
1733  END SELECT
1734  CASE(2)
1735  basis%DERIVATIVE_ORDER_INDEX_INV(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,1), &
1736  & basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,2),1,localnodeidx)=derivativeidx
1737  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,1))
1738  CASE(no_part_deriv)
1739  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,2))
1740  CASE(no_part_deriv)
1741  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=no_part_deriv
1742  CASE(first_part_deriv)
1743  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s2
1744  CASE DEFAULT
1745  CALL flagerror("Invalid derivative order index.",err,error,*999)
1746  END SELECT
1747  CASE(first_part_deriv)
1748  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,2))
1749  CASE(no_part_deriv)
1750  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s1
1751  CASE(first_part_deriv)
1752  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s1_s2
1753  CASE DEFAULT
1754  CALL flagerror("Invalid derivative order index.",err,error,*999)
1755  END SELECT
1756  CASE DEFAULT
1757  CALL flagerror("Invalid derivative order index.",err,error,*999)
1758  END SELECT
1759  CASE(3)
1760  basis%DERIVATIVE_ORDER_INDEX_INV(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,1), &
1761  & basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,2), &
1762  & basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,3),localnodeidx)=derivativeidx
1763  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,1))
1764  CASE(no_part_deriv)
1765  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,2))
1766  CASE(no_part_deriv)
1767  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,3))
1768  CASE(no_part_deriv)
1769  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=no_part_deriv
1770  CASE(first_part_deriv)
1771  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s3
1772  CASE DEFAULT
1773  CALL flagerror("Invalid derivative order index.",err,error,*999)
1774  END SELECT
1775  CASE(first_part_deriv)
1776  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,3))
1777  CASE(no_part_deriv)
1778  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s2
1779  CASE(first_part_deriv)
1780  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s2_s3
1781  CASE DEFAULT
1782  CALL flagerror("Invalid derivative order index.",err,error,*999)
1783  END SELECT
1784  CASE DEFAULT
1785  CALL flagerror("Invalid derivative order index.",err,error,*999)
1786  END SELECT
1787  CASE(first_part_deriv)
1788  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,2))
1789  CASE(no_part_deriv)
1790  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,3))
1791  CASE(no_part_deriv)
1792  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s1
1793  CASE(first_part_deriv)
1794  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s1_s3
1795  CASE DEFAULT
1796  CALL flagerror("Invalid derivative order index.",err,error,*999)
1797  END SELECT
1798  CASE(first_part_deriv)
1799  SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx,3))
1800  CASE(no_part_deriv)
1801  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s1_s2
1802  CASE(first_part_deriv)
1803  basis%PARTIAL_DERIVATIVE_INDEX(derivativeidx,localnodeidx)=part_deriv_s1_s2_s3
1804  CASE DEFAULT
1805  CALL flagerror("Invalid derivative order index.",err,error,*999)
1806  END SELECT
1807  CASE DEFAULT
1808  CALL flagerror("Invalid derivative order index.",err,error,*999)
1809  END SELECT
1810  CASE DEFAULT
1811  CALL flagerror("Invalid derivative order index.",err,error,*999)
1812  END SELECT
1813  CASE DEFAULT
1814  CALL flagerror("Invalid number of xi direcions.",err,error,*999)
1815  END SELECT
1816  ENDDO !derivativeIdx
1817  ENDDO !localNodeIdx
1818 
1819  !Set up the line information
1820  SELECT CASE(basis%NUMBER_OF_XI)
1821  CASE(1) !1 xi directions
1822  numberoflocallines=1
1823  basis%NUMBER_OF_LOCAL_LINES=1
1824  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(numberoflocallines),stat=err)
1825  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local line.",err,error,*999)
1826  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=basis%NUMBER_OF_NODES_XIC(1)
1827  ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(numberoflocallines),stat=err)
1828  IF(err/=0) CALL flagerror("Could not allocate local line xi direction.",err,error,*999)
1829  basis%LOCAL_LINE_XI_DIRECTION(1)=1
1830  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1),numberoflocallines),stat=err)
1831  IF(err/=0) CALL flagerror("Could not allocate node numbers in local line.",err,error,*999)
1832  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1),numberoflocallines),stat=err)
1833  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local line.",err,error,*999)
1834  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=no_part_deriv
1835  ALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1)**2,numberoflocallines),stat=err)
1836  IF(err/=0) CALL flagerror("Could not allocate element parameters in local line.",err,error,*999)
1837  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE=1
1838  locallineparameter=0
1839  DO localnodeidx2=1,basis%NUMBER_OF_NODES_XIC(1)
1840  DO localnodeidx1=1,basis%NUMBER_OF_NODES
1841  IF(basis%NODE_POSITION_INDEX(localnodeidx1,1)==localnodeidx2) THEN
1842  basis%NODE_NUMBERS_IN_LOCAL_LINE(localnodeidx2,1)=localnodeidx1
1843  DO derivativeidx=1,basis%NUMBER_OF_DERIVATIVES(localnodeidx2)
1844  locallineparameter=locallineparameter+1
1845  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(locallineparameter,1)=basis%ELEMENT_PARAMETER_INDEX( &
1846  & derivativeidx,localnodeidx1)
1847  IF(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnodeidx2,1)==first_part_deriv) THEN
1848  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(localnodeidx2,1)=derivativeidx
1849  EXIT
1850  ENDIF
1851  ENDDO !derivativeIdx
1852  EXIT
1853  ENDIF
1854  ENDDO !localNodeIdx2
1855  ENDDO !localNodeIdx1
1856  CASE(2) !2 xi directions
1857  !Determine the maximum node extent of the basis
1858  maximumnodeextent(1)=maxval(basis%NODE_POSITION_INDEX(:,1))
1859  maximumnodeextent(2)=maxval(basis%NODE_POSITION_INDEX(:,2))
1860  !Allocate and calculate the lines
1861  numberoflocallines=4-basis%NUMBER_OF_COLLAPSED_XI
1862  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(numberoflocallines),stat=err)
1863  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local line.",err,error,*999)
1864  basis%NUMBER_OF_NODES_IN_LOCAL_LINE=0
1865  ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(numberoflocallines),stat=err)
1866  IF(err/=0) CALL flagerror("Could not allocate local line xi direction.",err,error,*999)
1867  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),numberoflocallines),stat=err)
1868  IF(err/=0) CALL flagerror("Could not allocate node numbers in local line",err,error,*999)
1869  basis%NODE_NUMBERS_IN_LOCAL_LINE=0
1870  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),numberoflocallines),stat=err)
1871  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local line.",err,error,*999)
1872  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=no_part_deriv
1873  ALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC)**2,numberoflocallines),stat=err)
1874  IF(err/=0) CALL flagerror("Could not allocate element parameters in local line.",err,error,*999)
1875  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE=1
1876  ALLOCATE(basis%LOCAL_XI_NORMAL(numberoflocallines),stat=err)
1877  IF(err/=0) CALL flagerror("Could not allocate local xi normal.",err,error,*999)
1878  !Find the lines
1879  basis%NUMBER_OF_LOCAL_LINES=0
1880  DO xiidx1=1,2
1881  xiidx2=other_xi_directions2(xiidx1)
1882  !We are looking for lines in the xiIdx1 direction from the direction of xiIdx1=0
1883  !Loop over the element extremes in the xiIdx2 direction
1884  DO localnodeidx2=1,maximumnodeextent(xiidx2),maximumnodeextent(xiidx2)-1
1885  nodecount=0
1886  specialnodecount=0
1887  nodesinline=0
1888  DO localnodeidx1=1,basis%NUMBER_OF_NODES
1889  IF(basis%COLLAPSED_XI(xiidx1)/=basis_not_collapsed) THEN
1890  !The current xi direction, xiIdx1, is in a degenerate plane
1891  IF(basis%COLLAPSED_XI(xiidx2)==basis_xi_collapsed) THEN
1892  !The other xi direction is collapsed (must be the case)
1893  IF(basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi0) THEN !Collapsed at the xi=0 end
1894  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.OR. &
1895  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==1) THEN
1896  nodecount=nodecount+1
1897  nodesinline(nodecount)=localnodeidx1
1898  ENDIF
1899  ELSE !Collapsed at the xi=1 end
1900  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2) THEN
1901  nodecount=nodecount+1
1902  nodesinline(nodecount)=localnodeidx1
1903  ELSE IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==maximumnodeextent(xiidx1)) THEN
1904  IF(xiidx1<2) THEN !Special case - put the collapsed node at the end of the line
1905  specialnodecount=specialnodecount+1
1906  nodesinline(maximumnodeextent(xiidx1))=localnodeidx1
1907  ELSE
1908  nodecount=nodecount+1
1909  nodesinline(nodecount)=localnodeidx1
1910  ENDIF
1911  ENDIF
1912  ENDIF
1913  ELSE
1914  !The current xi direction must be collapsed
1915  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2) THEN
1916  nodecount=nodecount+1
1917  nodesinline(nodecount)=localnodeidx1
1918  ENDIF
1919  ENDIF
1920  ELSE
1921  !The current xi direction, xiIdx1, is not involved in any collapsed (degenerate) planes
1922  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2) THEN
1923  nodecount=nodecount+1
1924  nodesinline(nodecount)=localnodeidx1
1925  ENDIF
1926  ENDIF
1927  ENDDO !nn1
1928  IF((nodecount+specialnodecount)>1) THEN !More than one node so it is a proper line
1929  basis%NUMBER_OF_LOCAL_LINES=basis%NUMBER_OF_LOCAL_LINES+1
1930  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)=nodecount+specialnodecount
1931  basis%NODE_NUMBERS_IN_LOCAL_LINE(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES), &
1932  & basis%NUMBER_OF_LOCAL_LINES)=nodesinline(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES))
1933  locallineparameter=0
1934  DO locallinenodeidx=1,basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)
1935  DO derivativeidx=1,basis%NUMBER_OF_DERIVATIVES(basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
1936  & basis%NUMBER_OF_LOCAL_LINES))
1937  IF(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
1938  & basis%NUMBER_OF_LOCAL_LINES),xiidx2)==no_part_deriv) THEN
1939  locallineparameter=locallineparameter+1
1940  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(locallineparameter,basis%NUMBER_OF_LOCAL_LINES)= &
1941  & basis%ELEMENT_PARAMETER_INDEX(derivativeidx,basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
1942  & basis%NUMBER_OF_LOCAL_LINES))
1943  IF(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
1944  & basis%NUMBER_OF_LOCAL_LINES),xiidx1)==first_part_deriv) THEN
1945  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx,basis%NUMBER_OF_LOCAL_LINES)=derivativeidx
1946  ENDIF
1947  ENDIF
1948  ENDDO !derivativeIdx
1949  ENDDO !localLineNodeIdx
1950  basis%LOCAL_LINE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_LINES)=xiidx1
1951  IF(localnodeidx2==1) THEN
1952  basis%LOCAL_XI_NORMAL(basis%NUMBER_OF_LOCAL_LINES)=-xiidx2
1953  ELSE
1954  basis%LOCAL_XI_NORMAL(basis%NUMBER_OF_LOCAL_LINES)=xiidx2
1955  ENDIF
1956  ENDIF
1957  ENDDO !localNodeIdx2
1958  ENDDO !localNodeIdx1
1959  CASE(3) !3 xi directions
1960  !Determine the maximum node extent of the basis
1961  maximumnodeextent(1)=maxval(basis%NODE_POSITION_INDEX(:,1))
1962  maximumnodeextent(2)=maxval(basis%NODE_POSITION_INDEX(:,2))
1963  maximumnodeextent(3)=maxval(basis%NODE_POSITION_INDEX(:,3))
1964  !Allocate and calculate the lines
1965  IF(basis%NUMBER_OF_COLLAPSED_XI==1) THEN
1966  numberoflocallines=9
1967  numberoflocalfaces=5
1968  ELSE IF(basis%NUMBER_OF_COLLAPSED_XI==2) THEN
1969  numberoflocallines=8
1970  numberoflocalfaces=5
1971  ELSE
1972  numberoflocallines=12
1973  numberoflocalfaces=6
1974  ENDIF
1975  basis%NUMBER_OF_LOCAL_FACES=numberoflocalfaces
1976 
1977  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(numberoflocallines),stat=err)
1978  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local line.",err,error,*999)
1979  basis%NUMBER_OF_NODES_IN_LOCAL_LINE=0
1980 
1981  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_FACE(numberoflocalfaces),stat=err)
1982  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local face.",err,error,*999)
1983  basis%NUMBER_OF_NODES_IN_LOCAL_FACE=0
1984 
1985  ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(numberoflocallines),stat=err)
1986  IF(err/=0) CALL flagerror("Could not allocate local line xi direction.",err,error,*999)
1987 
1988  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),numberoflocallines),stat=err)
1989  IF(err/=0) CALL flagerror("Could not allocate node numbers in local line.",err,error,*999)
1990  basis%NODE_NUMBERS_IN_LOCAL_LINE=0
1991 
1992  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),numberoflocallines),stat=err)
1993  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local line.",err,error,*999)
1994  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=no_part_deriv
1995 
1996  ALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC)**2,numberoflocallines),stat=err)
1997  IF(err/=0) CALL flagerror("Could not allocate element parameters in local line.",err,error,*999)
1998  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE=1
1999 
2000  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:basis%MAXIMUM_NUMBER_OF_DERIVATIVES, &
2001  & maxval(basis%NUMBER_OF_NODES_XIC)**2,numberoflocalfaces),stat=err)
2002  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local face.",err,error,*999)
2003  basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE=no_part_deriv
2004  basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0,:,:)=1
2005 
2006  ALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_FACE(maxval(basis%NUMBER_OF_NODES_XIC)**2* &
2007  & basis%MAXIMUM_NUMBER_OF_DERIVATIVES,numberoflocalfaces),stat=err)
2008  IF(err/=0) CALL flagerror("Could not allocate element parameters in local face.",err,error,*999)
2009  basis%ELEMENT_PARAMETERS_IN_LOCAL_FACE=1
2010 
2011  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_FACE(max(maximumnodeextent(2)*maximumnodeextent(3), &
2012  & maximumnodeextent(3)*maximumnodeextent(1),maximumnodeextent(2)*maximumnodeextent(1)), &
2013  & numberoflocalfaces),stat=err)
2014  IF(err/=0) CALL flagerror("Could not allocate node numbers in local face.",err,error,*999)
2015  basis%NODE_NUMBERS_IN_LOCAL_FACE=0
2016 
2017  ALLOCATE(basis%LOCAL_XI_NORMAL(numberoflocalfaces),stat=err)
2018  IF(err/=0) CALL flagerror("Could not allocate local xi normal.",err,error,*999)
2019 
2020  ALLOCATE(basis%LOCAL_FACE_XI_DIRECTION(numberoflocalfaces),stat=err)
2021  IF(err/=0) CALL flagerror("Could not allocate local face xi direction.",err,error,*999)
2022 
2023  !Find the lines and faces
2024  basis%NUMBER_OF_LOCAL_LINES=0
2025  DO xiidx1=1,3
2026  xiidx2=other_xi_directions3(xiidx1,2,1)
2027  xiidx3=other_xi_directions3(xiidx1,3,1)
2028  !We are looking for lines going in the xiIdx1 direction, starting from xiIdx1=0.
2029  DO localnodeidx3=1,maximumnodeextent(xiidx3),maximumnodeextent(xiidx3)-1
2030  DO localnodeidx2=1,maximumnodeextent(xiidx2),maximumnodeextent(xiidx2)-1
2031  nodecount=0
2032  specialnodecount=0
2033  nodesinline=0
2034  !Iterate over nodes in the line of interest
2035  DO localnodeidx1=1,basis%NUMBER_OF_NODES
2036  IF(basis%COLLAPSED_XI(xiidx1)/=basis_not_collapsed) THEN
2037  !The current xi direction, xiIdx1, is involved in a collapsed (degenerate) plane
2038  IF(basis%COLLAPSED_XI(xiidx2)==basis_xi_collapsed.AND.basis%COLLAPSED_XI(xiidx3)==basis_xi_collapsed) THEN
2039  !Both of the other two xi directions are collapsed
2040  IF(basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi0) THEN !Collapsed at the xi=0 end
2041  IF((basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.OR. &
2042  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==1).AND. &
2043  & (basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3.OR. &
2044  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==1)) THEN
2045  nodecount=nodecount+1
2046  nodesinline(nodecount)=localnodeidx1
2047  ENDIF
2048  ELSE !Collapsed at the xi=1 end
2049  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2050  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2051  nodecount=nodecount+1
2052  nodesinline(nodecount)=localnodeidx1
2053  ELSE IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==maximumnodeextent(xiidx1)) THEN
2054  IF(xiidx1<3) THEN !Special case - put the collapsed node at the end of the line
2055  specialnodecount=specialnodecount+1
2056  nodesinline(maximumnodeextent(xiidx1))=localnodeidx1
2057  ELSE
2058  nodecount=nodecount+1
2059  nodesinline(nodecount)=localnodeidx1
2060  ENDIF
2061  ENDIF
2062  ENDIF
2063  ELSE
2064  IF(basis%COLLAPSED_XI(xiidx2)==basis_xi_collapsed) THEN
2065  !The other xiIdx2 xi direction is collapsed
2066  IF(basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi0) THEN !Collapsed at the xi=0 end
2067  IF((basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.OR. &
2068  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==1).AND. &
2069  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2070  nodecount=nodecount+1
2071  nodesinline(nodecount)=localnodeidx1
2072  ENDIF
2073  ELSE IF(basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi1) THEN !Collapsed at the xi=1 end
2074  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2075  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2076  nodecount=nodecount+1
2077  nodesinline(nodecount)=localnodeidx1
2078  ELSE IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==maximumnodeextent(xiidx1).AND. &
2079  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2080  IF(xiidx1<xiidx2) THEN !Special case - put the collapsed node at the end of the line
2081  specialnodecount=specialnodecount+1
2082  nodesinline(maximumnodeextent(xiidx1))=localnodeidx1
2083  ELSE
2084  nodecount=nodecount+1
2085  nodesinline(nodecount)=localnodeidx1
2086  ENDIF
2087  ENDIF
2088  ELSE
2089  !Not collapsed at a xi end
2090  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2091  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2092  nodecount=nodecount+1
2093  nodesinline(nodecount)=localnodeidx1
2094  ENDIF
2095  ENDIF
2096  ELSE IF(basis%COLLAPSED_XI(xiidx3)==basis_xi_collapsed) THEN
2097  !The other xiIdx3 xi direction is collapsed
2098  IF(basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi0) THEN !Collapsed at the xi=0 end
2099  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2100  & (basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3.OR. &
2101  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==1)) THEN
2102  nodecount=nodecount+1
2103  nodesinline(nodecount)=localnodeidx1
2104  ENDIF
2105  ELSE IF(basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi1) THEN !Collapsed at the xi=1 end
2106  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2107  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2108  nodecount=nodecount+1
2109  nodesinline(nodecount)=localnodeidx1
2110  ELSE IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx1)==maximumnodeextent(xiidx1).AND. &
2111  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2) THEN
2112  IF(xiidx1<xiidx3) THEN !Special case - put the collapsed node at the end of the line
2113  specialnodecount=specialnodecount+1
2114  nodesinline(maximumnodeextent(xiidx1))=localnodeidx1
2115  ELSE
2116  nodecount=nodecount+1
2117  nodesinline(nodecount)=localnodeidx1
2118  ENDIF
2119  ENDIF
2120  ELSE
2121  !Not collapsed at a xi end
2122  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2123  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2124  nodecount=nodecount+1
2125  nodesinline(nodecount)=localnodeidx1
2126  ENDIF
2127  ENDIF
2128  ELSE
2129  !The current xi must be collapsed
2130  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2131  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2132  nodecount=nodecount+1
2133  nodesinline(nodecount)=localnodeidx1
2134  ENDIF
2135  ENDIF
2136  ENDIF
2137  ELSE
2138  !The current xi direction, xiIdx1, is not involved in any collapsed (degenerate) planes
2139  IF(basis%NODE_POSITION_INDEX(localnodeidx1,xiidx2)==localnodeidx2.AND. &
2140  & basis%NODE_POSITION_INDEX(localnodeidx1,xiidx3)==localnodeidx3) THEN
2141  nodecount=nodecount+1
2142  nodesinline(nodecount)=localnodeidx1
2143  ENDIF
2144  ENDIF
2145  ENDDO !localNodeIdx1
2146  IF((nodecount+specialnodecount)>1) THEN !More than one node so it is a proper line
2147  basis%NUMBER_OF_LOCAL_LINES=basis%NUMBER_OF_LOCAL_LINES+1
2148  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)=nodecount+specialnodecount
2149  basis%NODE_NUMBERS_IN_LOCAL_LINE(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES), &
2150  & basis%NUMBER_OF_LOCAL_LINES)=nodesinline(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES))
2151  locallineparameter=0
2152  DO locallinenodeidx=1,basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)
2153  DO derivativeidx=1,basis%NUMBER_OF_DERIVATIVES(basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
2154  & basis%NUMBER_OF_LOCAL_LINES))
2155  IF(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
2156  & basis%NUMBER_OF_LOCAL_LINES),xiidx2)==no_part_deriv.AND. &
2157  & basis%DERIVATIVE_ORDER_INDEX(derivativeidx,basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
2158  & basis%NUMBER_OF_LOCAL_LINES),xiidx3)==no_part_deriv) THEN
2159  locallineparameter=locallineparameter+1
2160  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(locallineparameter,basis%NUMBER_OF_LOCAL_LINES)= &
2161  & basis%ELEMENT_PARAMETER_INDEX(derivativeidx,basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
2162  & basis%NUMBER_OF_LOCAL_LINES))
2163  IF(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,basis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx, &
2164  & basis%NUMBER_OF_LOCAL_LINES),xiidx1)==first_part_deriv) THEN
2165  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx,basis%NUMBER_OF_LOCAL_LINES)=derivativeidx
2166  ENDIF
2167  ENDIF
2168  ENDDO !derivativeIdx
2169  ENDDO !localLineNodeIdx
2170  basis%LOCAL_LINE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_LINES)=xiidx1
2171  ENDIF
2172  ENDDO !localNodeIdx2
2173  ENDDO !localNodeIdx3
2174  ENDDO !xiIdx1
2175 
2176  !Find the local nodes and derivatives in each face and the local face xi direction
2177  localfaceidx=0
2178  !Loop over the -'ve and +'ve xi direction
2179  DO directionidx=-1,1,2
2180  !Loop over the three xi directions
2181  DO xiidx1=1,3
2182  !xiIdx1 is the +/- face normal direction. xiIdx2 and xiIdx3 are the xi directions in the face.
2183  xiidx2=other_xi_directions3(xiidx1,2,1)
2184  xiidx3=other_xi_directions3(xiidx1,3,1)
2185 
2186  IF(directionidx==1) THEN
2187  !The +'ve xi direction
2188  localnodeidx1=maximumnodeextent(xiidx1)
2189  !Compute if the face in the +xiIdx1 direction is collapsed.
2190  collapsedface=basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi1
2191  ELSE
2192  !The -'ve xi direction
2193  localnodeidx1=1
2194  !Compute if the face in the +xiIdx1 direction is collapsed.
2195  collapsedface=basis%COLLAPSED_XI(xiidx1)==basis_collapsed_at_xi0
2196  ENDIF
2197  localnodecount=0
2198  IF(.NOT.collapsedface) THEN
2199  !If the face has not been collapsed
2200  localfaceidx=localfaceidx+1
2201  !Loop over the local nodes in the face
2202  DO localnodeidx3=1,maximumnodeextent(xiidx2)
2203  DO localnodeidx2=1,maximumnodeextent(xiidx3)
2204  IF(xiidx1==1) THEN
2205  localnodeidx=basis%NODE_POSITION_INDEX_INV(localnodeidx1,localnodeidx2,localnodeidx3,1)
2206  ELSE IF(xiidx1==2) THEN
2207  localnodeidx=basis%NODE_POSITION_INDEX_INV(localnodeidx2,localnodeidx1,localnodeidx3,1)
2208  ELSE
2209  localnodeidx=basis%NODE_POSITION_INDEX_INV(localnodeidx2,localnodeidx3,localnodeidx1,1)
2210  ENDIF
2211  IF(all(basis%NODE_NUMBERS_IN_LOCAL_FACE(1:localnodecount,localfaceidx)/=localnodeidx)) THEN
2212  !The node hasn't been collapsed
2213  localnodecount=localnodecount+1
2214  basis%NODE_NUMBERS_IN_LOCAL_FACE(localnodecount,localfaceidx)=localnodeidx
2215  ENDIF
2216  ENDDO !localNodeIdx3
2217  ENDDO !localNodexIdx2
2218  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(localfaceidx)=localnodecount
2219  basis%LOCAL_FACE_XI_DIRECTION(localfaceidx)=directionidx*xiidx1
2220  !Compute derivatives and element parameters in the face
2221  localfaceparameter=0
2222  DO localnodeidx=1,basis%NUMBER_OF_NODES_IN_LOCAL_FACE(localfaceidx)
2223  localnode=basis%NODE_NUMBERS_IN_LOCAL_FACE(localnodeidx,localfaceidx)
2224  localfacederivative=0
2225  DO derivativeidx=1,basis%NUMBER_OF_DERIVATIVES(localnode)
2226  IF(basis%DERIVATIVE_ORDER_INDEX(derivativeidx,localnode,xiidx1)==no_part_deriv) THEN
2227  localfaceparameter=localfaceparameter+1
2228  localfacederivative=localfacederivative+1
2229  basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(localfacederivative,localnodeidx,localfaceidx)=derivativeidx
2230  basis%ELEMENT_PARAMETERS_IN_LOCAL_FACE(localfaceparameter,localfaceidx)= &
2231  & basis%ELEMENT_PARAMETER_INDEX(derivativeidx,localnode)
2232  ENDIF
2233  ENDDO !derivativeIdx
2234  basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0,localnodeidx,localfaceidx)=localfacederivative
2235  ENDDO !localNodeIdx
2236  ENDIF
2237  ENDDO !xiIdx1
2238  ENDDO !directionIdx
2239  CASE DEFAULT
2240  CALL flagerror("Invalid number of xi directions.",err,error,*999)
2241  END SELECT
2242 
2243  CALL basis_quadrature_create(basis,err,error,*999)
2244 
2245  ELSE
2246  CALL flagerror("Basis is not a Lagrange Hermite tensor product basis.",err,error,*999)
2247  ENDIF
2248  ELSE
2249  CALL flagerror("Basis is not associated.",err,error,*999)
2250  ENDIF
2251 
2252  exits("Basis_LHTPBasisCreate")
2253  RETURN
2254 999 IF(ALLOCATED(nodeatcollapse)) DEALLOCATE(nodeatcollapse)
2255  errorsexits("Basis_LHTPBasisCreate",err,error)
2256  RETURN 1
2257  END SUBROUTINE basis_lhtpbasiscreate
2258 
2259  !
2260  !
2261  !================================================================================================================================
2262  !
2263 
2265  FUNCTION basis_lhtp_basis_evaluate_dp(BASIS,NODE_NUMBER,DERIVATIVE_NUMBER,PARTIAL_DERIV_INDEX,XI,ERR,ERROR)
2266 
2267  !Argument variables
2268  TYPE(basis_type), POINTER :: BASIS
2269  INTEGER(INTG), INTENT(IN) :: NODE_NUMBER
2270  INTEGER(INTG), INTENT(IN) :: DERIVATIVE_NUMBER
2271  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIV_INDEX
2272  REAL(DP), INTENT(IN) :: XI(:)
2273  INTEGER(INTG), INTENT(OUT) :: ERR
2274  TYPE(varying_string), INTENT(OUT) :: ERROR
2275  !Function variable
2276  REAL(DP) :: BASIS_LHTP_BASIS_EVALUATE_DP
2277  !Local variables
2278  INTEGER(INTG) :: ni,nn
2279  REAL(DP) :: SUM
2280  TYPE(varying_string) :: LOCAL_ERROR
2281 
2282  enters("BASIS_LHTP_BASIS_EVALUATE_DP",err,error,*999)
2283 
2284  basis_lhtp_basis_evaluate_dp=1.0_dp
2285  IF(ASSOCIATED(basis)) THEN
2286  DO ni=1,basis%NUMBER_OF_XI
2287  IF(basis%NODE_AT_COLLAPSE(node_number).AND.basis%COLLAPSED_XI(ni)==basis_xi_collapsed) THEN
2288  !We are at a collapsed node in the collapsed xi direction. Sum the basis functions in the collapsed xi direction.
2289  sum=0.0_dp
2290  SELECT CASE(basis%INTERPOLATION_TYPE(ni))
2291  CASE(basis_lagrange_interpolation)
2292  SELECT CASE(basis%INTERPOLATION_ORDER(ni))
2293  CASE(basis_linear_interpolation_order)
2294  DO nn=1,2
2295  sum=sum+lagrange_linear_evaluate(nn,partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2296  ENDDO !nn
2297  CASE(basis_quadratic_interpolation_order)
2298  DO nn=1,3
2299  sum=sum+lagrange_quadratic_evaluate(nn,partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2300  ENDDO !nn
2301  CASE(basis_cubic_interpolation_order)
2302  DO nn=1,4
2303  sum=sum+lagrange_cubic_evaluate(nn,partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2304  ENDDO !nn
2305  CASE DEFAULT
2306  local_error="Interpolation order value "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(ni),"*",err,error))// &
2307  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid"
2308  CALL flagerror(local_error,err,error,*999)
2309  END SELECT
2310  CASE(basis_hermite_interpolation)
2311  SELECT CASE(basis%INTERPOLATION_ORDER(ni))
2312  CASE(basis_cubic_interpolation_order)
2313  DO nn=1,2
2314  sum=sum+hermite_cubic_evaluate(nn,basis%DERIVATIVE_ORDER_INDEX(derivative_number,node_number,ni), &
2315  & partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2316  ENDDO !nn
2317  CASE DEFAULT
2318  local_error="Interpolation order value "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(ni),"*",err,error))// &
2319  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid"
2320  CALL flagerror(local_error,err,error,*999)
2321  END SELECT
2322  IF(err/=0) GOTO 999
2323  CASE DEFAULT
2324  local_error="Interpolation type value "//trim(number_to_vstring(basis%INTERPOLATION_TYPE(ni),"*",err,error))// &
2325  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid"
2326  CALL flagerror(local_error,err,error,*999)
2327  END SELECT
2328  basis_lhtp_basis_evaluate_dp=basis_lhtp_basis_evaluate_dp*sum
2329  ELSE
2330  SELECT CASE(basis%INTERPOLATION_TYPE(ni))
2331  CASE(basis_lagrange_interpolation)
2332  SELECT CASE(basis%INTERPOLATION_ORDER(ni))
2333  CASE(basis_linear_interpolation_order)
2334  basis_lhtp_basis_evaluate_dp=basis_lhtp_basis_evaluate_dp* &
2335  & lagrange_linear_evaluate(basis%NODE_POSITION_INDEX(node_number,ni), &
2336  & partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2337  CASE(basis_quadratic_interpolation_order)
2338  basis_lhtp_basis_evaluate_dp=basis_lhtp_basis_evaluate_dp* &
2339  & lagrange_quadratic_evaluate(basis%NODE_POSITION_INDEX(node_number,ni), &
2340  & partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2341  CASE(basis_cubic_interpolation_order)
2342  basis_lhtp_basis_evaluate_dp=basis_lhtp_basis_evaluate_dp* &
2343  & lagrange_cubic_evaluate(basis%NODE_POSITION_INDEX(node_number,ni), &
2344  & partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2345  CASE DEFAULT
2346  local_error="Interpolation order value "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(ni),"*",err,error))// &
2347  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid"
2348  CALL flagerror(local_error,err,error,*999)
2349  END SELECT
2350  CASE(basis_hermite_interpolation)
2351  SELECT CASE(basis%INTERPOLATION_ORDER(ni))
2352  CASE(basis_quadratic1_interpolation_order)
2353  basis_lhtp_basis_evaluate_dp=basis_lhtp_basis_evaluate_dp* &
2354  & hermite_quadratic_evaluate(basis%NODE_POSITION_INDEX(node_number,ni), &
2355  & basis%DERIVATIVE_ORDER_INDEX(derivative_number,node_number,ni), &
2356  & partial_derivative_index(partial_deriv_index,ni),1,xi(ni),err,error)
2357  CASE(basis_quadratic2_interpolation_order)
2358  basis_lhtp_basis_evaluate_dp=basis_lhtp_basis_evaluate_dp* &
2359  & hermite_quadratic_evaluate(basis%NODE_POSITION_INDEX(node_number,ni), &
2360  & basis%DERIVATIVE_ORDER_INDEX(derivative_number,node_number,ni), &
2361  & partial_derivative_index(partial_deriv_index,ni),2,xi(ni),err,error)
2362  CASE(basis_cubic_interpolation_order)
2363  basis_lhtp_basis_evaluate_dp=basis_lhtp_basis_evaluate_dp* &
2364  & hermite_cubic_evaluate(basis%NODE_POSITION_INDEX(node_number,ni), &
2365  & basis%DERIVATIVE_ORDER_INDEX(derivative_number,node_number,ni), &
2366  & partial_derivative_index(partial_deriv_index,ni),xi(ni),err,error)
2367  CASE DEFAULT
2368  local_error="Interpolation order value "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(ni),"*",err,error))// &
2369  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid"
2370  CALL flagerror(local_error,err,error,*999)
2371  END SELECT
2372  IF(err/=0) GOTO 999
2373  CASE DEFAULT
2374  local_error="Interpolation type value "//trim(number_to_vstring(basis%INTERPOLATION_TYPE(ni),"*",err,error))// &
2375  & " for xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid"
2376  CALL flagerror(local_error,err,error,*999)
2377  END SELECT
2378  ENDIF
2379  ENDDO !ni
2380  ELSE
2381  CALL flagerror("Basis is not associated",err,error,*999)
2382  ENDIF
2383 
2384  exits("BASIS_LHTP_BASIS_EVALUATE_DP")
2385  RETURN
2386 999 errorsexits("BASIS_LHTP_BASIS_EVALUATE_DP",err,error)
2387  RETURN
2388  END FUNCTION basis_lhtp_basis_evaluate_dp
2389 
2390  !
2391  !================================================================================================================================
2392  !
2393 
2396  SUBROUTINE basis_lhtp_family_create(BASIS,ERR,ERROR,*)
2397 
2398  !Argument variables
2399  TYPE(basis_type), POINTER :: BASIS
2400  INTEGER(INTG), INTENT(OUT) :: ERR
2401  TYPE(varying_string), INTENT(OUT) :: ERROR
2402  !Local Variables
2403  INTEGER(INTG) :: DUMMY_ERR,ni,ni2,FACE_XI(2)
2404  LOGICAL :: LINE_BASIS_DONE,FACE_BASIS_DONE
2405  TYPE(basis_type), POINTER :: NEW_SUB_BASIS
2406  TYPE(varying_string) :: DUMMY_ERROR
2407 
2408  NULLIFY(new_sub_basis)
2409 
2410  enters("BASIS_LHTP_FAMILY_CREATE",err,error,*999)
2411 
2412  IF(ASSOCIATED(basis)) THEN
2413  !Create the main (parent) basis
2414  CALL basis_lhtpbasiscreate(basis,err,error,*999)
2415  IF(basis%NUMBER_OF_XI>1) THEN
2416  !Create the line bases as sub-basis types
2417  ALLOCATE(basis%LINE_BASES(basis%NUMBER_OF_XI),stat=err)
2418  IF(err/=0) CALL flagerror("Could not allocate basis line bases",err,error,*999)
2419  DO ni=1,basis%NUMBER_OF_XI
2420  line_basis_done=.false.
2421  NULLIFY(new_sub_basis)
2422  DO ni2=1,ni-1
2423  IF(basis%INTERPOLATION_XI(ni2)==basis%INTERPOLATION_XI(ni).AND. &
2424  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni2)==basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)) THEN
2425  line_basis_done=.true.
2426  EXIT
2427  ENDIF
2428  ENDDO !ni2
2429  IF(line_basis_done) THEN
2430  basis%LINE_BASES(ni)%PTR=>basis%LINE_BASES(ni2)%PTR
2431  ELSE
2432  !Create the new sub-basis
2433  CALL basis_sub_basis_create(basis,1,[ni],new_sub_basis,err,error,*999)
2434  !Fill in the basis information
2435  CALL basis_lhtpbasiscreate(new_sub_basis,err,error,*999)
2436  basis%LINE_BASES(ni)%PTR=>new_sub_basis
2437  ENDIF
2438  ENDDO !ni
2439  IF(basis%NUMBER_OF_XI>2) THEN
2440  !Set up face basis functions
2441  ALLOCATE(basis%FACE_BASES(basis%NUMBER_OF_XI),stat=err)
2442  IF(err/=0) CALL flagerror("Could not allocate basis face bases",err,error,*999)
2443  DO ni=1,basis%NUMBER_OF_XI
2444  !Determine the face xi directions that lie in this xi direction
2445  face_xi(1)=other_xi_directions3(ni,2,1)
2446  face_xi(2)=other_xi_directions3(ni,3,1)
2447  face_basis_done=.false.
2448  NULLIFY(new_sub_basis)
2449  DO ni2=1,ni-1
2450 ! FACE_XI2(1)=OTHER_XI_DIRECTIONS3(ni2,2,1)
2451 ! FACE_XI2(2)=OTHER_XI_DIRECTIONS3(ni2,3,1)
2452 
2453 !!TODO FIX THIS
2454  !Going to disable the test below, as it results in error in collapsed elements and doesn't save much time
2455 ! IF(BASIS%INTERPOLATION_XI(FACE_XI2(1))==BASIS%INTERPOLATION_XI(FACE_XI(1)).AND. &
2456 ! & BASIS%INTERPOLATION_XI(FACE_XI2(2))==BASIS%INTERPOLATION_XI(FACE_XI(2)).AND. &
2457 ! & BASIS%QUADRATURE%NUMBER_OF_GAUSS_XI(FACE_XI2(1))==BASIS%QUADRATURE%NUMBER_OF_GAUSS_XI(FACE_XI(1)).AND. &
2458 ! & BASIS%QUADRATURE%NUMBER_OF_GAUSS_XI(FACE_XI2(2))==BASIS%QUADRATURE%NUMBER_OF_GAUSS_XI(FACE_XI(1))) THEN
2459 ! FACE_BASIS_DONE=.TRUE.
2460 ! EXIT
2461  ! ENDIF
2462 
2463  ENDDO !ni2
2464  IF(face_basis_done) THEN
2465  basis%FACE_BASES(ni)%PTR=>basis%FACE_BASES(ni2)%PTR
2466  ELSE
2467  !Create the new sub-basis
2468  CALL basis_sub_basis_create(basis,2,[face_xi(1),face_xi(2)],new_sub_basis,err,error,*999)
2469  !Fill in the basis information
2470  CALL basis_lhtpbasiscreate(new_sub_basis,err,error,*999)
2471  new_sub_basis%LINE_BASES(1)%PTR=>basis%LINE_BASES(face_xi(1))%PTR
2472  new_sub_basis%LINE_BASES(2)%PTR=>basis%LINE_BASES(face_xi(2))%PTR
2473  basis%FACE_BASES(ni)%PTR=>new_sub_basis
2474  ENDIF
2475  ENDDO !ni
2476  ELSE
2477  ALLOCATE(basis%FACE_BASES(1),stat=err)
2478  IF(err/=0) CALL flagerror("Could not allocate basis face bases",err,error,*999)
2479  basis%FACE_BASES(1)%PTR=>basis
2480  ENDIF
2481  ELSE
2482  ALLOCATE(basis%LINE_BASES(1),stat=err)
2483  IF(err/=0) CALL flagerror("Could not allocate basis line bases",err,error,*999)
2484  basis%LINE_BASES(1)%PTR=>basis
2485  NULLIFY(basis%FACE_BASES)
2486  ENDIF
2487  ELSE
2488  CALL flagerror("Basis is not associated",err,error,*999)
2489  ENDIF
2490 
2491  exits("BASIS_LHTP_FAMILY_CREATE")
2492  RETURN
2493 999 IF(ASSOCIATED(new_sub_basis)) CALL basis_family_destroy(new_sub_basis%USER_NUMBER,new_sub_basis%FAMILY_NUMBER, &
2494  & dummy_err,dummy_error,*998)
2495 998 errorsexits("BASIS_LHTP_FAMILY_CREATE",err,error)
2496  RETURN 1
2497  END SUBROUTINE basis_lhtp_family_create
2498 
2499  !
2500  !================================================================================================================================
2501  !
2502 
2505  SUBROUTINE basis_radial_family_create(BASIS,ERR,ERROR,*)
2506 
2507  !Argument variables
2508  TYPE(basis_type), POINTER :: BASIS
2509  INTEGER(INTG), INTENT(OUT) :: ERR
2510  TYPE(varying_string), INTENT(OUT) :: ERROR
2511  !Local Variables
2512  TYPE(varying_string) :: DUMMY_ERROR
2513 
2514  enters("BASIS_RADIAL_FAMILY_CREATE",err,error,*999)
2515 
2516  IF(ASSOCIATED(basis)) THEN
2517  !Create the main (parent) basis
2518  CALL flagerror("Not implemented.",err,error,*999)
2519  ELSE
2520  CALL flagerror("Basis is not associated",err,error,*999)
2521  ENDIF
2522 
2523  exits("BASIS_RADIAL_FAMILY_CREATE")
2524  RETURN
2525 999 errorsexits("BASIS_RADIAL_FAMILY_CREATE",err,error)
2526  RETURN 1
2527 
2528  END SUBROUTINE basis_radial_family_create
2529 
2530  !
2531  !================================================================================================================================
2532  !
2533 
2535  SUBROUTINE basis_local_node_xi_calculate(BASIS,LOCAL_NODE_NUMBER,XI,ERR,ERROR,*)
2536 
2537  !Argument variables
2538  TYPE(basis_type), POINTER :: BASIS
2539  INTEGER(INTG), INTENT(IN) :: LOCAL_NODE_NUMBER
2540  REAL(DP), INTENT(OUT) :: XI(:)
2541  INTEGER(INTG), INTENT(OUT) :: ERR
2542  TYPE(varying_string), INTENT(OUT) :: ERROR
2543  !Local Variables
2544  INTEGER(INTG) :: xi_idx
2545  TYPE(varying_string) :: LOCAL_ERROR
2546 
2547  enters("BASIS_LOCAL_NODE_XI_CALCULATE",err,error,*999)
2548 
2549  IF(ASSOCIATED(basis)) THEN
2550  IF(basis%BASIS_FINISHED) THEN
2551  IF(local_node_number>0.AND.local_node_number<=basis%NUMBER_OF_NODES) THEN
2552  IF(SIZE(xi,1)>=basis%NUMBER_OF_XI) THEN
2553  SELECT CASE(basis%TYPE)
2554  CASE(basis_lagrange_hermite_tp_type)
2555  DO xi_idx=1,basis%NUMBER_OF_XI
2556  xi(xi_idx)=REAL(basis%node_position_index(local_node_number,xi_idx)-1,dp)/ &
2557  & REAL(BASIS%NUMBER_OF_NODES_XIC(xi_idx)-1,DP)
2558  ENDDO !xi_idx
2559  CASE(basis_simplex_type)
2560  DO xi_idx=1,basis%NUMBER_OF_XI
2561  xi(xi_idx)=REAL(BASIS%NUMBER_OF_NODES_XIC(xi_idx)-BASIS%NODE_POSITION_INDEX(LOCAL_NODE_NUMBER,xi_idx),DP)/ &
2562  & REAL(BASIS%NUMBER_OF_NODES_XIC(xi_idx)-1,DP)
2563  ENDDO !xi_idx
2564  CASE(basis_serendipity_type)
2565  CALL flagerror("Not implemented.",err,error,*999)
2566  CASE(basis_auxilliary_type)
2567  CALL flagerror("Not implemented.",err,error,*999)
2568  CASE(basis_b_spline_tp_type)
2569  CALL flagerror("Not implemented.",err,error,*999)
2570  CASE(basis_fourier_lagrange_hermite_tp_type)
2571  CALL flagerror("Not implemented.",err,error,*999)
2572  CASE(basis_extended_lagrange_tp_type)
2573  CALL flagerror("Not implemented.",err,error,*999)
2574  CASE DEFAULT
2575  local_error="The basis type of "//trim(number_to_vstring(basis%TYPE,"*",err,error))// &
2576  & " is invalid."
2577  CALL flagerror(local_error,err,error,*999)
2578  END SELECT
2579  ELSE
2580  local_error="The size of the specified xic array of "//trim(number_to_vstring(SIZE(xi,1),"*",err,error))// &
2581  & " is invalid. The size of the xi array must be >= "// &
2582  & trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))//"."
2583  CALL flagerror(local_error,err,error,*999)
2584  ENDIF
2585  ELSE
2586  local_error="The specified local node number of "//trim(number_to_vstring(local_node_number,"*",err,error))// &
2587  & " is invalid. The local node number must be > 0 and <= "// &
2588  & trim(number_to_vstring(basis%NUMBER_OF_NODES,"*",err,error))//"."
2589  CALL flagerror(local_error,err,error,*999)
2590  ENDIF
2591  ELSE
2592  CALL flagerror("Basis has not been finished.",err,error,*999)
2593  ENDIF
2594  ELSE
2595  CALL flagerror("Basis is not associated",err,error,*999)
2596  ENDIF
2597 
2598  exits("BASIS_LOCAL_NODE_XI_CALCULATE")
2599  RETURN
2600 999 errorsexits("BASIS_LOCAL_NODE_XI_CALCULATE",err,error)
2601  RETURN 1
2602  END SUBROUTINE basis_local_node_xi_calculate
2603 
2604  !
2605  !================================================================================================================================
2606  !
2607 
2609  SUBROUTINE basis_number_of_local_nodes_get(BASIS,NUMBER_OF_LOCAL_NODES,ERR,ERROR,*)
2610 
2611  !Argument variables
2612  TYPE(basis_type), POINTER :: BASIS
2613  INTEGER(INTG), INTENT(OUT) :: NUMBER_OF_LOCAL_NODES
2614  INTEGER(INTG), INTENT(OUT) :: ERR
2615  TYPE(varying_string), INTENT(OUT) :: ERROR
2616  !Local Variables
2617 
2618  enters("BASIS_NUMBER_OF_LOCAL_NODES_GET",err,error,*999)
2619 
2620  IF(ASSOCIATED(basis)) THEN
2621  number_of_local_nodes=basis%NUMBER_OF_NODES
2622  ELSE
2623  CALL flagerror("Basis is not associated",err,error,*999)
2624  ENDIF
2625 
2626  exits("BASIS_NUMBER_OF_LOCAL_NODES_GET")
2627  RETURN
2628 999 errorsexits("BASIS_NUMBER_OF_LOCAL_NODES_GET",err,error)
2629  RETURN 1
2630  END SUBROUTINE basis_number_of_local_nodes_get
2631 
2632  !
2633  !================================================================================================================================
2634  !
2635 
2637  SUBROUTINE basis_number_of_xi_get(BASIS,NUMBER_OF_XI,ERR,ERROR,*)
2638 
2639  !Argument variables
2640  TYPE(basis_type), POINTER :: BASIS
2641  INTEGER(INTG), INTENT(OUT) :: NUMBER_OF_XI
2642  INTEGER(INTG), INTENT(OUT) :: ERR
2643  TYPE(varying_string), INTENT(OUT) :: ERROR
2644  !Local Variables
2645 
2646  enters("BASIS_NUMBER_OF_XI_GET",err,error,*999)
2647 
2648  IF(ASSOCIATED(basis)) THEN
2649  IF(basis%BASIS_FINISHED) THEN
2650  number_of_xi=basis%NUMBER_OF_XI
2651  ELSE
2652  CALL flagerror("Basis has not been finished.",err,error,*999)
2653  ENDIF
2654  ELSE
2655  CALL flagerror("Basis is not associated.",err,error,*999)
2656  ENDIF
2657 
2658  exits("BASIS_NUMBER_OF_XI_GET")
2659  RETURN
2660 999 errorsexits("BASIS_NUMBER_OF_XI_GET",err,error)
2661  RETURN
2662  END SUBROUTINE basis_number_of_xi_get
2663 
2664  !
2665  !================================================================================================================================
2666  !
2667 
2669  SUBROUTINE basis_number_of_xi_set_number(USER_NUMBER,NUMBER_OF_XI,ERR,ERROR,*)
2670 
2671  !Argument variables
2672  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
2673  INTEGER(INTG), INTENT(IN) :: NUMBER_OF_XI
2674  INTEGER(INTG), INTENT(OUT) :: ERR
2675  TYPE(varying_string), INTENT(OUT) :: ERROR
2676  !Local Variables
2677  TYPE(basis_type), POINTER :: BASIS
2678 
2679  enters("BASIS_NUMBER_OF_XI_SET_NUMBER",err,error,*999)
2680 
2681  CALL basis_user_number_find(user_number,basis,err,error,*999)
2682  CALL basis_number_of_xi_set(basis,number_of_xi,err,error,*999)
2683 
2684  exits("BASIS_NUMBER_OF_XI_SET_NUMBER")
2685  RETURN
2686 999 errorsexits("BASIS_NUMBER_OF_XI_SET_NUMBER",err,error)
2687  RETURN 1
2688  END SUBROUTINE basis_number_of_xi_set_number
2689 
2690  !
2691  !================================================================================================================================
2692  !
2693 
2695  SUBROUTINE basis_number_of_xi_set_ptr(BASIS,NUMBER_OF_XI,ERR,ERROR,*)
2696 
2697  !Argument variables
2698  TYPE(basis_type), POINTER :: BASIS
2699  INTEGER(INTG), INTENT(IN) :: NUMBER_OF_XI
2700  INTEGER(INTG), INTENT(OUT) :: ERR
2701  TYPE(varying_string), INTENT(OUT) :: ERROR
2702  !Local Variables
2703  INTEGER(INTG) :: OLD_INTERPOLATION_XI(3),OLD_NUMBER_OF_GAUSS_XI(3),OLD_COLLAPSED_XI(3)
2704  TYPE(varying_string) :: LOCAL_ERROR
2705 
2706  enters("BASIS_NUMBER_OF_XI_SET_PTR",err,error,*999)
2707 
2708  IF(ASSOCIATED(basis)) THEN
2709  IF(basis%BASIS_FINISHED) THEN
2710  CALL flagerror("Basis has been finished",err,error,*999)
2711  ELSE
2712  SELECT CASE(basis%TYPE)
2713  CASE(basis_lagrange_hermite_tp_type)
2714  IF(number_of_xi>0.AND.number_of_xi<4) THEN
2715  IF(basis%NUMBER_OF_XI/=number_of_xi) THEN
2716  !Reallocate the basis information arrays that depend on the number of xi directions
2717  old_interpolation_xi=basis%INTERPOLATION_XI
2718  old_collapsed_xi=basis%COLLAPSED_XI
2719  DEALLOCATE(basis%INTERPOLATION_XI)
2720  DEALLOCATE(basis%COLLAPSED_XI)
2721  ALLOCATE(basis%INTERPOLATION_XI(number_of_xi),stat=err)
2722  IF(err/=0) CALL flagerror("Could not allocate interpolation type",err,error,*999)
2723  ALLOCATE(basis%COLLAPSED_XI(number_of_xi),stat=err)
2724  IF(err/=0) CALL flagerror("Could not allocate collapsed xi",err,error,*999)
2725  IF(number_of_xi>basis%NUMBER_OF_XI) THEN
2726  basis%INTERPOLATION_XI(1:basis%NUMBER_OF_XI)=old_interpolation_xi(1:basis%NUMBER_OF_XI)
2727  basis%INTERPOLATION_XI(basis%NUMBER_OF_XI+1:number_of_xi)=old_interpolation_xi(1)
2728  basis%COLLAPSED_XI(1:basis%NUMBER_OF_XI)=old_collapsed_xi(1:basis%NUMBER_OF_XI)
2729  basis%COLLAPSED_XI(basis%NUMBER_OF_XI+1:number_of_xi)=old_collapsed_xi(1)
2730  ELSE
2731  basis%INTERPOLATION_XI(1:number_of_xi)=old_interpolation_xi(1:number_of_xi)
2732  basis%COLLAPSED_XI(1:number_of_xi)=old_collapsed_xi(1:number_of_xi)
2733  ENDIF
2734 
2735  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
2736  old_number_of_gauss_xi=basis%QUADRATURE%NUMBER_OF_GAUSS_XI
2737  DEALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI)
2738  ALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(number_of_xi),stat=err)
2739  IF(err/=0) CALL flagerror("Could not allocate number of Gauss xi",err,error,*999)
2740  IF(number_of_xi>basis%NUMBER_OF_XI) THEN
2741  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI)=old_number_of_gauss_xi(1:basis%NUMBER_OF_XI)
2742  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(basis%NUMBER_OF_XI+1:number_of_xi)=old_number_of_gauss_xi(1)
2743  ELSE
2744  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:number_of_xi)=old_number_of_gauss_xi(1:number_of_xi)
2745  ENDIF
2746  ENDIF
2747  basis%NUMBER_OF_XI=number_of_xi
2748  ENDIF
2749  ELSE
2750  local_error="Invalid number of xi directions specified ("// &
2751  & trim(number_to_vstring(number_of_xi,"*",err,error))// &
2752  & ") for a Lagrange-Hermite basis. You must specify between 1 and 3 xi directions"
2753  CALL flagerror(local_error,err,error,*999)
2754  ENDIF
2755  CASE(basis_simplex_type)
2756  IF(number_of_xi>1.AND.number_of_xi<4) THEN
2757  IF(basis%NUMBER_OF_XI/=number_of_xi) THEN
2758  !Reallocate the basis information arrays that depend on the number of xi directions
2759  old_interpolation_xi=basis%INTERPOLATION_XI
2760  old_collapsed_xi=basis%COLLAPSED_XI
2761  DEALLOCATE(basis%INTERPOLATION_XI)
2762  DEALLOCATE(basis%COLLAPSED_XI)
2763  ALLOCATE(basis%INTERPOLATION_XI(number_of_xi),stat=err)
2764  IF(err/=0) CALL flagerror("Could not allocate interpolation type",err,error,*999)
2765  ALLOCATE(basis%COLLAPSED_XI(number_of_xi),stat=err)
2766  IF(err/=0) CALL flagerror("Could not allocate collapsed xi.",err,error,*999)
2767  IF(number_of_xi>basis%NUMBER_OF_XI) THEN
2768  basis%INTERPOLATION_XI(1:basis%NUMBER_OF_XI)=old_interpolation_xi(1:basis%NUMBER_OF_XI)
2769  basis%INTERPOLATION_XI(basis%NUMBER_OF_XI+1:number_of_xi)=old_interpolation_xi(1)
2770  basis%COLLAPSED_XI(1:basis%NUMBER_OF_XI)=old_collapsed_xi(1:basis%NUMBER_OF_XI)
2771  basis%COLLAPSED_XI(basis%NUMBER_OF_XI+1:number_of_xi)=old_collapsed_xi(1)
2772  ELSE
2773  basis%INTERPOLATION_XI(1:number_of_xi)=old_interpolation_xi(1:number_of_xi)
2774  basis%COLLAPSED_XI(1:number_of_xi)=old_collapsed_xi(1:number_of_xi)
2775  ENDIF
2776  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
2777  old_number_of_gauss_xi=basis%QUADRATURE%NUMBER_OF_GAUSS_XI
2778  DEALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI)
2779  ALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(number_of_xi),stat=err)
2780  IF(err/=0) CALL flagerror("Could not allocate number of Gauss xi",err,error,*999)
2781  IF(number_of_xi>basis%NUMBER_OF_XI) THEN
2782  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI)=old_number_of_gauss_xi(1:basis%NUMBER_OF_XI)
2783  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(basis%NUMBER_OF_XI+1:number_of_xi)=old_number_of_gauss_xi(1)
2784  ELSE
2785  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:number_of_xi)=old_number_of_gauss_xi(1:number_of_xi)
2786  ENDIF
2787  ENDIF
2788  basis%NUMBER_OF_XI=number_of_xi
2789  ENDIF
2790  ELSE
2791  local_error="Invalid number of xi directions specified ("// &
2792  & trim(number_to_vstring(number_of_xi,"*",err,error))// &
2793  & ") for a simplex basis. You must specify between 2 and 3 xi directions"
2794  CALL flagerror(local_error,err,error,*999)
2795  ENDIF
2796  CASE DEFAULT
2797  CALL flagerror("Basis type invalid or not implemented",err,error,*999)
2798  END SELECT
2799  ENDIF
2800  ELSE
2801  CALL flagerror("Basis is not associated",err,error,*999)
2802  ENDIF
2803 
2804  exits("BASIS_NUMBER_OF_XI_SET_PTR")
2805  RETURN
2806 999 errorsexits("BASIS_NUMBER_OF_XI_SET_PTR",err,error)
2807  RETURN 1
2808  END SUBROUTINE basis_number_of_xi_set_ptr
2809 
2810  !
2811  !================================================================================================================================
2812  !
2813 
2815  SUBROUTINE basis_quadrature_create(BASIS,ERR,ERROR,*)
2816 
2817  !Argument variables
2818  TYPE(basis_type), POINTER :: BASIS
2819  INTEGER(INTG), INTENT(OUT) :: ERR
2820  TYPE(varying_string), INTENT(OUT) :: ERROR
2821  !Local Variables
2822  INTEGER(INTG) :: scheme_idx,i,j,k,MAX_NUM_GAUSS,ng,ni,nk,nn,ns,nu,NUM_GAUSS_1,NUM_GAUSS_2,NUM_GAUSS_3
2823  REAL(DP) :: XI(3),GSX(4,20),GSW(20)
2824  REAL(DP), ALLOCATABLE :: POSITIONS(:,:),POSITIONS_MATRIX(:,:,:,:),WEIGHTS(:,:)
2825  TYPE(quadrature_scheme_type), POINTER :: NEW_SCHEME,SCHEME
2826  TYPE(quadrature_scheme_ptr_type), POINTER :: NEW_SCHEMES(:)
2827  TYPE(varying_string) :: LOCAL_ERROR
2828  INTEGER(INTG) :: MAX_NUM_FACE_GAUSS,face_idx,NORMAL,FACE_XI(2),numberOfFaceXiCoordinates
2829 
2830  NULLIFY(new_scheme)
2831  NULLIFY(new_schemes)
2832 
2833  enters("BASIS_QUADRATURE_CREATE",err,error,*999)
2834 
2835  IF(ASSOCIATED(basis)) THEN
2836  IF(ASSOCIATED(basis%QUADRATURE%SCHEMES)) THEN
2837  local_error="The quadrature schemes on basis number "//trim(number_to_vstring(basis%USER_NUMBER,"*",err,error))// &
2838  & " are already associated"
2839  CALL flagerror(local_error,err,error,*998)
2840  ELSE
2841 !!TODO: \todo Sort this properly by having a create values cache.
2842  !Reset the basis quadrature -
2843  !CALL BASIS_QUADRATURE_FINALISE(BASIS,ERR,ERROR,*999) ! Kumar - I think this is not correct as it
2844  !Initialise the basis quadrature ! resets the quadrature scheme already set.
2845  !CALL BASIS_QUADRATURE_INITIALISE(BASIS,ERR,ERROR,*999) !
2846  SELECT CASE(basis%QUADRATURE%TYPE)
2847  CASE(basis_gauss_legendre_quadrature)
2848  !Allocate one scheme and add it to the list of schemes
2849  ALLOCATE(new_scheme,stat=err)
2850  IF(err/=0) CALL flagerror("Could not allocate new quadrature scheme",err,error,*999)
2851  new_scheme%QUADRATURE=>basis%QUADRATURE
2852  basis%QUADRATURE%NUMBER_OF_SCHEMES=1
2853  ALLOCATE(new_schemes(basis%QUADRATURE%NUMBER_OF_SCHEMES),stat=err)
2854  IF(err/=0) CALL flagerror("Could not allocate new quadratures scheme",err,error,*999)
2855  new_schemes(1)%PTR=>new_scheme
2856  basis%QUADRATURE%SCHEMES=>new_schemes
2857  !Set up the quadrature scheme map
2858  ALLOCATE(basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_number_of_quadrature_scheme_types),stat=err)
2859  IF(err/=0) CALL flagerror("Could not allocate quadrature scheme map",err,error,*999)
2860  DO scheme_idx=1,basis_number_of_quadrature_scheme_types
2861  NULLIFY(basis%QUADRATURE%QUADRATURE_SCHEME_MAP(scheme_idx)%PTR)
2862  ENDDO !scheme_idx
2863  basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR=>new_scheme
2864  !Set up the gauss point arrays
2865  new_scheme%NUMBER_OF_GAUSS=1
2866  max_num_gauss=-1
2867  DO ni=1,basis%NUMBER_OF_XI
2868  new_scheme%NUMBER_OF_GAUSS=new_scheme%NUMBER_OF_GAUSS*basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)
2869  IF(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)>max_num_gauss) max_num_gauss=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)
2870  ENDDO !ni
2871  ALLOCATE(new_scheme%GAUSS_POSITIONS(basis%NUMBER_OF_XI_COORDINATES,new_scheme%NUMBER_OF_GAUSS),stat=err)
2872  IF(err/=0) CALL flagerror("Could not allocate Gauss positions",err,error,*999)
2873  ALLOCATE(new_scheme%GAUSS_WEIGHTS(new_scheme%NUMBER_OF_GAUSS),stat=err)
2874  IF(err/=0) CALL flagerror("Could not allocate Gauss weights",err,error,*999)
2875  ALLOCATE(new_scheme%GAUSS_BASIS_FNS(basis%NUMBER_OF_ELEMENT_PARAMETERS,basis%NUMBER_OF_PARTIAL_DERIVATIVES, &
2876  & new_scheme%NUMBER_OF_GAUSS),stat=err)
2877  IF(err/=0) CALL flagerror("Could not allocate Gauss basis functions",err,error,*999)
2878  ALLOCATE(weights(max_num_gauss,3),stat=err)
2879  IF(err/=0) CALL flagerror("Could not allocate weights",err,error,*999)
2880  ALLOCATE(positions(max_num_gauss,3),stat=err)
2881  IF(err/=0) CALL flagerror("Could not allocate positions",err,error,*999)
2882  ALLOCATE(positions_matrix(max_num_gauss,max_num_gauss,max_num_gauss,3),stat=err)
2883  IF(err/=0) CALL flagerror("Could not allocate positions matrix",err,error,*999)
2884  weights=1.0_dp
2885  positions=0.0_dp
2886  positions_matrix=0.0_dp
2887  DO ni=1,basis%NUMBER_OF_XI
2888  CALL gauss_legendre(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),0.0_dp,1.0_dp, &
2889  & positions(1:basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),ni), &
2890  & weights(1:basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),ni),err,error,*999)
2891  ENDDO !ni
2892  SELECT CASE(basis%NUMBER_OF_XI)
2893  CASE(1)
2894  num_gauss_1=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1)
2895  num_gauss_2=1
2896  num_gauss_3=1
2897  CASE(2)
2898  num_gauss_1=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1)
2899  num_gauss_2=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(2)
2900  num_gauss_3=1
2901  CASE(3)
2902  num_gauss_1=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1)
2903  num_gauss_2=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(2)
2904  num_gauss_3=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(3)
2905  CASE DEFAULT
2906  CALL flagerror("Invalid number of xi directions",err,error,*999)
2907  END SELECT
2908  DO k=1,num_gauss_3
2909  DO j=1,num_gauss_2
2910  DO i=1,num_gauss_1
2911  positions_matrix(i,j,k,1)=positions(i,1)
2912  positions_matrix(i,j,k,2)=positions(j,2)
2913  positions_matrix(i,j,k,3)=positions(k,3)
2914  xi(1:basis%NUMBER_OF_XI)=positions_matrix(i,j,k,1:basis%NUMBER_OF_XI)
2915  ng=i+(j-1+(k-1)*num_gauss_2)*num_gauss_1
2916  new_scheme%GAUSS_WEIGHTS(ng)=weights(i,1)*weights(j,2)*weights(k,3)
2917  new_scheme%GAUSS_POSITIONS(1:basis%NUMBER_OF_XI_COORDINATES,ng)=xi(1:basis%NUMBER_OF_XI_COORDINATES)
2918  ns=0
2919  DO nn=1,basis%NUMBER_OF_NODES
2920  DO nk=1,basis%NUMBER_OF_DERIVATIVES(nn)
2921  ns=ns+1
2922  DO nu=1,basis%NUMBER_OF_PARTIAL_DERIVATIVES
2923  SELECT CASE(basis%TYPE)
2924  CASE(basis_lagrange_hermite_tp_type)
2925  new_scheme%GAUSS_BASIS_FNS(ns,nu,ng)=basis_lhtp_basis_evaluate(basis,nn,nk,nu,xi,err,error)
2926  IF(err/=0) GOTO 999
2927  CASE DEFAULT
2928  CALL flagerror("Not implemented",err,error,*999)
2929  END SELECT
2930  ENDDO !nu
2931  ENDDO !nk
2932  ENDDO !nn
2933  ENDDO !i
2934  ENDDO !j
2935  ENDDO !k
2936  !Create face quadrature scheme, if requested
2937  IF(basis%QUADRATURE%EVALUATE_FACE_GAUSS) THEN
2938  IF(basis%NUMBER_OF_XI==3) THEN
2939  !Find maximum number of face gauss points and allocate the arrays
2940  max_num_face_gauss=product(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI))
2941  max_num_face_gauss=max_num_face_gauss/minval(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI))
2942  ALLOCATE(new_scheme%NUMBER_OF_FACE_GAUSS(basis%NUMBER_OF_LOCAL_FACES),stat=err)
2943  IF(err/=0) CALL flagerror("Could not allocate number of face gauss",err,error,*999)
2944  ALLOCATE(new_scheme%FACE_GAUSS_POSITIONS(basis%NUMBER_OF_XI_COORDINATES,max_num_face_gauss, &
2945  & basis%NUMBER_OF_LOCAL_FACES),stat=err)
2946  IF(err/=0) CALL flagerror("Could not allocate face Gauss positions",err,error,*999)
2947  ALLOCATE(new_scheme%FACE_GAUSS_WEIGHTS(max_num_face_gauss,basis%NUMBER_OF_LOCAL_FACES),stat=err)
2948  IF(err/=0) CALL flagerror("Could not allocate face Gauss weights",err,error,*999)
2949  ALLOCATE(new_scheme%FACE_GAUSS_BASIS_FNS(basis%NUMBER_OF_ELEMENT_PARAMETERS,basis%NUMBER_OF_PARTIAL_DERIVATIVES, &
2950  & max_num_face_gauss,basis%NUMBER_OF_LOCAL_FACES),stat=err)
2951  IF(err/=0) CALL flagerror("Could not allocate face Gauss basis function values array",err,error,*999)
2952  !Zero them out just to be safe
2953  new_scheme%FACE_GAUSS_POSITIONS=0.0_dp
2954  new_scheme%FACE_GAUSS_WEIGHTS=0.0_dp
2955  new_scheme%FACE_GAUSS_BASIS_FNS=0.0_dp
2956  !Populate face_gauss_positions, weights, basis_fn
2957  DO face_idx=1,basis%NUMBER_OF_LOCAL_FACES
2958  !What's the normal?
2959  normal=basis%LOCAL_FACE_XI_DIRECTION(face_idx)
2960  IF(normal<0_intg) THEN
2961  xi(abs(normal))=0.0_dp
2962  ELSE
2963  xi(abs(normal))=1.0_dp
2964  ENDIF
2965  normal=abs(normal)
2966  face_xi=[other_xi_directions3(normal,2,1), other_xi_directions3(normal,3,1)]
2967  !How many gauss points are in this face?
2968  new_scheme%NUMBER_OF_FACE_GAUSS(face_idx)=product(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(face_xi))
2969  ng=0_intg
2970  DO j=1,basis%QUADRATURE%NUMBER_OF_GAUSS_XI(face_xi(2))
2971  xi(face_xi(2))=positions(j,face_xi(2))
2972  DO i=1,basis%QUADRATURE%NUMBER_OF_GAUSS_XI(face_xi(1))
2973  xi(face_xi(1))=positions(i,face_xi(1))
2974  ng=ng+1_intg
2975  !Gauss point xi and weights first
2976  new_scheme%FACE_GAUSS_WEIGHTS(ng,face_idx)=weights(i,face_xi(1))*weights(j,face_xi(2))
2977  new_scheme%FACE_GAUSS_POSITIONS(1:3,ng,face_idx)=xi(1:3)
2978  !Evaluate basis fn values at the Gauss points now
2979  ns=0
2980  DO nn=1,basis%NUMBER_OF_NODES
2981  DO nk=1,basis%NUMBER_OF_DERIVATIVES(nn)
2982  ns=ns+1
2983  DO nu=1,basis%NUMBER_OF_PARTIAL_DERIVATIVES
2984  SELECT CASE(basis%TYPE)
2985  CASE(basis_lagrange_hermite_tp_type)
2986  new_scheme%FACE_GAUSS_BASIS_FNS(ns,nu,ng,face_idx)= &
2987  & basis_lhtp_basis_evaluate(basis,nn,nk,nu,xi,err,error)
2988  IF(err/=0) GOTO 999
2989  CASE DEFAULT
2990  CALL flagerror("Not implemented",err,error,*999)
2991  END SELECT
2992  ENDDO !nu
2993  ENDDO !nk
2994  ENDDO !nn
2995 
2996  ENDDO !i
2997  ENDDO !j
2998  ENDDO !face_idx
2999  ELSE
3000  CALL flagerror("Cannot evaluate face quadrature schemes for a non three dimensional element.",err,error,*999)
3001  ENDIF
3002  ENDIF
3003  !Clean up
3004  DEALLOCATE(weights)
3005  DEALLOCATE(positions)
3006  DEALLOCATE(positions_matrix)
3007  CASE(basis_gauss_laguerre_quadrature)
3008  CALL flagerror("Gauss Laguerre quadrature type not implemented.",err,error,*999)
3009  CASE(basis_guass_hermite_quadrature)
3010  CALL flagerror("Gauss Hermite quadrature type not implemented.",err,error,*999)
3011  CASE(basis_gauss_simplex_quadrature)
3012  !Allocate one scheme and add it to the list of schemes
3013  ALLOCATE(new_scheme,stat=err)
3014  IF(err/=0) CALL flagerror("Could not allocate new quadrature scheme.",err,error,*999)
3015  new_scheme%QUADRATURE=>basis%QUADRATURE
3016  basis%QUADRATURE%NUMBER_OF_SCHEMES=1
3017  ALLOCATE(new_schemes(basis%QUADRATURE%NUMBER_OF_SCHEMES),stat=err)
3018  IF(err/=0) CALL flagerror("Could not allocate new quadratures scheme.",err,error,*999)
3019  new_schemes(1)%PTR=>new_scheme
3020  basis%QUADRATURE%SCHEMES=>new_schemes
3021  !Set up the quadrature scheme map
3022  ALLOCATE(basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_number_of_quadrature_scheme_types),stat=err)
3023  IF(err/=0) CALL flagerror("Could not allocate quadrature scheme map.",err,error,*999)
3024  DO scheme_idx=1,basis_number_of_quadrature_scheme_types
3025  NULLIFY(basis%QUADRATURE%QUADRATURE_SCHEME_MAP(scheme_idx)%PTR)
3026  ENDDO !scheme_idx
3027  basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR=>new_scheme
3028  !Set up the gauss point arrays
3029  CALL gauss_simplex(basis%QUADRATURE%GAUSS_ORDER,basis%NUMBER_OF_XI_COORDINATES,new_scheme%NUMBER_OF_GAUSS,gsx,gsw, &
3030  & err,error,*999)
3031  ALLOCATE(new_scheme%GAUSS_POSITIONS(basis%NUMBER_OF_XI_COORDINATES,new_scheme%NUMBER_OF_GAUSS),stat=err)
3032  IF(err/=0) CALL flagerror("Could not allocate Gauss positions.",err,error,*999)
3033  ALLOCATE(new_scheme%GAUSS_WEIGHTS(new_scheme%NUMBER_OF_GAUSS),stat=err)
3034  IF(err/=0) CALL flagerror("Could not allocate Gauss weights.",err,error,*999)
3035  ALLOCATE(new_scheme%GAUSS_BASIS_FNS(basis%NUMBER_OF_ELEMENT_PARAMETERS,basis%NUMBER_OF_PARTIAL_DERIVATIVES, &
3036  & new_scheme%NUMBER_OF_GAUSS),stat=err)
3037  IF(err/=0) CALL flagerror("Could not allocate Gauss basis functions.",err,error,*999)
3038  new_scheme%GAUSS_POSITIONS(1:basis%NUMBER_OF_XI_COORDINATES,1:new_scheme%NUMBER_OF_GAUSS)= &
3039  & gsx(1:basis%NUMBER_OF_XI_COORDINATES,1:new_scheme%NUMBER_OF_GAUSS)
3040  new_scheme%GAUSS_WEIGHTS(1:new_scheme%NUMBER_OF_GAUSS)=gsw(1:new_scheme%NUMBER_OF_GAUSS)
3041  DO ng=1,new_scheme%NUMBER_OF_GAUSS
3042  ns=0
3043  DO nn=1,basis%NUMBER_OF_NODES
3044  DO nk=1,basis%NUMBER_OF_DERIVATIVES(nn)
3045  ns=ns+1
3046  DO nu=1,basis%NUMBER_OF_PARTIAL_DERIVATIVES
3047  SELECT CASE(basis%TYPE)
3048  CASE(basis_simplex_type)
3049  !Gauss positions are in area coordinates so call the simplex basis evaluate directly
3050  new_scheme%GAUSS_BASIS_FNS(ns,nu,ng)= &
3051  & basis_simplex_basis_evaluate(basis,nn,nu,new_scheme%GAUSS_POSITIONS(1:basis%NUMBER_OF_XI_COORDINATES,ng), &
3052  & err,error)
3053  IF(err/=0) GOTO 999
3054  CASE DEFAULT
3055  CALL flagerror("Not implemented.",err,error,*999)
3056  END SELECT
3057  ENDDO !nu
3058  ENDDO !nk
3059  ENDDO !nn
3060  ENDDO !ng
3061  !Create face quadrature scheme, if requested
3062  IF(basis%QUADRATURE%EVALUATE_FACE_GAUSS) THEN
3063  IF(basis%NUMBER_OF_XI==3) THEN
3064  !Find maximum number of face gauss points and allocate the arrays
3065  max_num_face_gauss=product(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI))
3066  max_num_face_gauss=max_num_face_gauss/minval(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI))
3067  ALLOCATE(new_scheme%NUMBER_OF_FACE_GAUSS(basis%NUMBER_OF_LOCAL_FACES),stat=err)
3068  IF(err/=0) CALL flagerror("Could not allocate number of face gauss",err,error,*999)
3069  ALLOCATE(new_scheme%FACE_GAUSS_POSITIONS(basis%NUMBER_OF_XI_COORDINATES,max_num_face_gauss, &
3070  & basis%NUMBER_OF_LOCAL_FACES),stat=err)
3071  IF(err/=0) CALL flagerror("Could not allocate face Gauss positions",err,error,*999)
3072  ALLOCATE(new_scheme%FACE_GAUSS_WEIGHTS(max_num_face_gauss,basis%NUMBER_OF_LOCAL_FACES),stat=err)
3073  IF(err/=0) CALL flagerror("Could not allocate face Gauss weights",err,error,*999)
3074  ALLOCATE(new_scheme%FACE_GAUSS_BASIS_FNS(basis%NUMBER_OF_ELEMENT_PARAMETERS,basis%NUMBER_OF_PARTIAL_DERIVATIVES, &
3075  & max_num_face_gauss,basis%NUMBER_OF_LOCAL_FACES),stat=err)
3076  IF(err/=0) CALL flagerror("Could not allocate face Gauss basis function values array",err,error,*999)
3077  !Zero them out just to be safe
3078  new_scheme%FACE_GAUSS_POSITIONS=0.0_dp
3079  new_scheme%FACE_GAUSS_WEIGHTS=0.0_dp
3080  new_scheme%FACE_GAUSS_BASIS_FNS=0.0_dp
3081  !Populate face_gauss_positions, weights, basis_fn
3082  DO face_idx=1,basis%NUMBER_OF_LOCAL_FACES
3083  !The number of face xi coordinates will be 3 for triangular face on a tet
3084  numberoffacexicoordinates = basis%NUMBER_OF_XI
3085  !Set up the gauss point arrays for the face
3086  CALL gauss_simplex(basis%QUADRATURE%GAUSS_ORDER,numberoffacexicoordinates, &
3087  & new_scheme%NUMBER_OF_FACE_GAUSS(face_idx),gsx,gsw,err,error,*999)
3088  IF(err/=0) CALL flagerror("Could not allocate Gauss basis functions",err,error,*999)
3089  new_scheme%FACE_GAUSS_POSITIONS(1:numberoffacexicoordinates,1:new_scheme%NUMBER_OF_FACE_GAUSS(face_idx), &
3090  & face_idx)=gsx(1:numberoffacexicoordinates,1:new_scheme%NUMBER_OF_FACE_GAUSS(face_idx))
3091  new_scheme%FACE_GAUSS_WEIGHTS(1:new_scheme%NUMBER_OF_FACE_GAUSS(face_idx),face_idx)= &
3092  & gsw(1:new_scheme%NUMBER_OF_FACE_GAUSS(face_idx))
3093 
3094  DO ng=1,new_scheme%NUMBER_OF_FACE_GAUSS(face_idx)
3095  ns=0
3096  DO nn=1,basis%NUMBER_OF_NODES_IN_LOCAL_FACE(face_idx)
3097  DO nk=1,basis%NUMBER_OF_DERIVATIVES(nn)
3098  ns=ns+1
3099  DO nu=1,basis%NUMBER_OF_PARTIAL_DERIVATIVES
3100  SELECT CASE(basis%TYPE)
3101  CASE(basis_simplex_type)
3102  new_scheme%FACE_GAUSS_BASIS_FNS(ns,nu,ng,face_idx)= &
3103  & basis_simplex_basis_evaluate(basis,nn,nu, &
3104  & new_scheme%FACE_GAUSS_POSITIONS(1:numberoffacexicoordinates,ng,face_idx),err,error)
3105  IF(err/=0) GOTO 999
3106  CASE DEFAULT
3107  CALL flagerror("Not implemented",err,error,*999)
3108  END SELECT
3109  ENDDO !nu
3110  ENDDO !nk
3111  ENDDO !nn
3112  ENDDO !ng
3113 
3114  ENDDO !face_idx
3115  ELSE
3116  CALL flagerror("Cannot evaluate face quadrature schemes for a non three dimensional element.",err,error,*999)
3117  ENDIF
3118  ENDIF
3119  CASE DEFAULT
3120  local_error="Quadrature type "//trim(number_to_vstring(basis%QUADRATURE%TYPE,"*",err,error))//" is invalid."
3121  CALL flagerror(local_error,err,error,*999)
3122  END SELECT
3123  ENDIF
3124  ELSE
3125  CALL flagerror("Basis is not associated.",err,error,*998)
3126  ENDIF
3127 
3128  IF(diagnostics1) THEN
3129  CALL write_string_value(diagnostic_output_type,"Quadrature type = ",basis%QUADRATURE%TYPE,err,error,*999)
3130  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_XI,3,3,basis%QUADRATURE%NUMBER_OF_GAUSS_XI, &
3131  & '(" Number of gauss points(ni):",3(X,I2))','(22X,3(X,I2))',err,error,*999)
3132  CALL write_string_value(diagnostic_output_type," Number of quadrature schemes = ",basis%QUADRATURE%NUMBER_OF_SCHEMES, &
3133  & err,error,*999)
3134  DO scheme_idx=1,basis%QUADRATURE%NUMBER_OF_SCHEMES
3135  scheme=>basis%QUADRATURE%SCHEMES(scheme_idx)%PTR
3136  CALL write_string_value(diagnostic_output_type," Scheme = ",scheme_idx,err,error,*999)
3137  CALL write_string_value(diagnostic_output_type," Total number of gauss points = ",scheme%NUMBER_OF_GAUSS, &
3138  & err,error,*999)
3139  IF(diagnostics2) THEN
3140  CALL write_string(diagnostic_output_type," Gauss point positions and weights:",err,error,*999)
3141  DO ng=1,scheme%NUMBER_OF_GAUSS
3142  CALL write_string_value(diagnostic_output_type," Gauss point = ",ng,err,error,*999)
3143  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_XI_COORDINATES,3,3,scheme%GAUSS_POSITIONS(:,ng), &
3144  & '(" position(ni) :",3(X,F12.4))','(26X,3(X,F12.4))',err,error,*999)
3145  CALL write_string_fmt_value(diagnostic_output_type," WEIGHT : ",scheme%GAUSS_WEIGHTS(ng), &
3146  & "(F12.4)",err,error,*999)
3147  ENDDO !ng
3148  ENDIF
3149  IF(diagnostics3) THEN
3150  CALL write_string(diagnostic_output_type," Basis functions evaluated at Gauss points:",err,error,*999)
3151  DO ng=1,scheme%NUMBER_OF_GAUSS
3152  CALL write_string_value(diagnostic_output_type," Gauss point = ",ng,err,error,*999)
3153  DO nu=1,basis%NUMBER_OF_PARTIAL_DERIVATIVES
3154  CALL write_string_value(diagnostic_output_type," Partial derivative number = ",nu,err,error,*999)
3155  CALL write_string_vector(diagnostic_output_type,1,1,basis%NUMBER_OF_ELEMENT_PARAMETERS,4,4, &
3156  & scheme%GAUSS_BASIS_FNS(:,nu,ng),'(" BASIS FNS(ns) :",4(X,F12.4))','(26X,4(X,F12.4))',err,error,*999)
3157  ENDDO !nu
3158  ENDDO !ng
3159  ENDIF
3160  ENDDO !scheme_idx
3161  ENDIF
3162 
3163  exits("BASIS_QUADRATURE_CREATE")
3164  RETURN
3165 999 IF(ASSOCIATED(new_scheme)) THEN
3166  IF(ALLOCATED(new_scheme%GAUSS_POSITIONS)) DEALLOCATE(new_scheme%GAUSS_POSITIONS)
3167  IF(ALLOCATED(new_scheme%GAUSS_WEIGHTS)) DEALLOCATE(new_scheme%GAUSS_WEIGHTS)
3168  IF(ALLOCATED(new_scheme%GAUSS_BASIS_FNS)) DEALLOCATE(new_scheme%GAUSS_BASIS_FNS)
3169  DEALLOCATE(new_scheme)
3170  ENDIF
3171  IF(ALLOCATED(weights)) DEALLOCATE(weights)
3172  IF(ALLOCATED(positions)) DEALLOCATE(positions)
3173  IF(ALLOCATED(positions_matrix)) DEALLOCATE(positions_matrix)
3174  IF(ASSOCIATED(new_schemes)) DEALLOCATE(new_schemes)
3175  NULLIFY(basis%QUADRATURE%SCHEMES)
3176 998 errorsexits("BASIS_QUADRATURE_CREATE",err,error)
3177  RETURN 1
3178 
3179  END SUBROUTINE basis_quadrature_create
3180 
3181  !
3182  !================================================================================================================================
3183  !
3184 
3186  SUBROUTINE basis_quadrature_destroy(QUADRATURE,ERR,ERROR,*)
3187 
3188  !Argument variables
3189  TYPE(quadrature_type), POINTER :: QUADRATURE
3190  INTEGER(INTG), INTENT(OUT) :: ERR
3191  TYPE(varying_string), INTENT(OUT) :: ERROR
3192  !Local Variables
3193 
3194  enters("BASIS_QUADRATURE_DESTROY",err,error,*999)
3195 
3196  IF(ASSOCIATED(quadrature)) THEN
3197  CALL flagerror("Not implemented.",err,error,*999)
3198  ELSE
3199  CALL flagerror("Basis is not associated",err,error,*999)
3200  ENDIF
3201 
3202  exits("BASIS_QUADRATURE_DESTROY")
3203  RETURN
3204 999 errorsexits("BASIS_QUADRATURE_DESTROY",err,error)
3205  RETURN 1
3206  END SUBROUTINE basis_quadrature_destroy
3207 
3208  !
3209  !================================================================================================================================
3210  !
3211 
3213  SUBROUTINE basis_quadrature_finalise(BASIS,ERR,ERROR,*)
3214 
3215  !Argument variables
3216  TYPE(basis_type), POINTER :: BASIS
3217  INTEGER(INTG), INTENT(OUT) :: ERR
3218  TYPE(varying_string), INTENT(OUT) :: ERROR
3219  !Local Variables
3220  INTEGER(INTG) :: scheme_idx
3221  TYPE(quadrature_scheme_type), POINTER :: SCHEME
3222 
3223  enters("BASIS_QUADRATURE_FINALISE",err,error,*999)
3224 
3225  IF(ASSOCIATED(basis)) THEN
3226  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3227  DO scheme_idx=1,basis%QUADRATURE%NUMBER_OF_SCHEMES
3228  scheme=>basis%QUADRATURE%SCHEMES(scheme_idx)%PTR
3229  !Destroy all scheme components
3230  IF (ASSOCIATED(scheme)) THEN
3231  IF(ALLOCATED(scheme%GAUSS_POSITIONS)) DEALLOCATE(scheme%GAUSS_POSITIONS)
3232  IF(ALLOCATED(scheme%GAUSS_WEIGHTS)) DEALLOCATE(scheme%GAUSS_WEIGHTS)
3233  IF(ALLOCATED(scheme%GAUSS_BASIS_FNS)) DEALLOCATE(scheme%GAUSS_BASIS_FNS)
3234  DEALLOCATE(scheme)
3235  ENDIF
3236  ENDDO !scheme_idx
3237  IF(ASSOCIATED(basis%QUADRATURE%SCHEMES)) DEALLOCATE(basis%QUADRATURE%SCHEMES)
3238  basis%QUADRATURE%NUMBER_OF_SCHEMES=0
3239  IF(ALLOCATED(basis%QUADRATURE%NUMBER_OF_GAUSS_XI)) DEALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI)
3240  NULLIFY(basis%QUADRATURE%BASIS)
3241  basis%QUADRATURE%TYPE=-1
3242  IF(ALLOCATED(basis%QUADRATURE%QUADRATURE_SCHEME_MAP)) DEALLOCATE(basis%QUADRATURE%QUADRATURE_SCHEME_MAP)
3243  ELSE
3244  CALL flagerror("Basis quadrature basis is not associated",err,error,*999)
3245  ENDIF
3246  ELSE
3247  CALL flagerror("Basis is not associated",err,error,*999)
3248  ENDIF
3249 
3250  exits("BASIS_QUADRATURE_FINALISE")
3251  RETURN
3252 999 errorsexits("BASIS_QUADRATURE_FINALISE",err,error)
3253  RETURN 1
3254  END SUBROUTINE basis_quadrature_finalise
3255 
3256  !
3257  !================================================================================================================================
3258  !
3259 
3261  SUBROUTINE basis_quadrature_initialise(BASIS,ERR,ERROR,*)
3262 
3263  !Argument variables
3264  TYPE(basis_type), POINTER :: BASIS
3265  INTEGER(INTG), INTENT(OUT) :: ERR
3266  TYPE(varying_string), INTENT(OUT) :: ERROR
3267  !Local Variables
3268  INTEGER(INTG) :: ni
3269  TYPE(varying_string) :: LOCAL_ERROR
3270 
3271  enters("BASIS_QUADRATURE_INITIALISE",err,error,*999)
3272 
3273  IF(ASSOCIATED(basis)) THEN
3274  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3275  local_error="Basis number "//trim(number_to_vstring(basis%USER_NUMBER,"*",err,error))// &
3276  & " already has a quadrature associated"
3277  CALL flagerror(local_error,err,error,*998)
3278  ELSE
3279  SELECT CASE(basis%TYPE)
3280  CASE(basis_lagrange_hermite_tp_type)
3281  !Set up a default Gauss Legendre quadrature
3282  basis%QUADRATURE%NUMBER_OF_SCHEMES=0
3283  NULLIFY(basis%QUADRATURE%SCHEMES)
3284  basis%QUADRATURE%BASIS=>basis
3285  basis%QUADRATURE%TYPE=basis_gauss_legendre_quadrature
3286  !Set up a default number of Gauss points appropriate for the given interpolation order in each direction.
3287  ALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(basis%NUMBER_OF_XI),stat=err)
3288  IF(err/=0) CALL flagerror("Could not allocate number of Gauss in each xi direction",err,error,*999)
3289  DO ni=1,basis%NUMBER_OF_XI
3290  SELECT CASE(basis%INTERPOLATION_XI(ni))
3291  CASE(basis_linear_lagrange_interpolation)
3292  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)=2
3293  CASE(basis_quadratic_lagrange_interpolation)
3294  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)=3
3295  CASE(basis_cubic_lagrange_interpolation)
3296  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)=4
3297  CASE(basis_quadratic1_hermite_interpolation,basis_quadratic2_hermite_interpolation)
3298  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)=3
3299  CASE(basis_cubic_hermite_interpolation)
3300  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)=4
3301  CASE DEFAULT
3302  local_error="Interpolation xi value "//trim(number_to_vstring(basis%INTERPOLATION_XI(ni),"*",err,error))// &
3303  & " in xi direction "//trim(number_to_vstring(ni,"*",err,error))//" is invalid"
3304  CALL flagerror(local_error,err,error,*999)
3305  END SELECT
3306  ENDDO !ni
3307  CASE(basis_simplex_type)
3308  !Set up a default quadrature
3309  basis%QUADRATURE%NUMBER_OF_SCHEMES=0
3310  NULLIFY(basis%QUADRATURE%SCHEMES)
3311  basis%QUADRATURE%BASIS=>basis
3312  basis%QUADRATURE%TYPE=basis_gauss_simplex_quadrature
3313  !Set up a default order appropriate for the given interpolation.
3314  ALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(basis%NUMBER_OF_XI),stat=err)
3315  IF(err/=0) CALL flagerror("Could not allocate number of Gauss in each xi direction",err,error,*999)
3316 !!TODO: \todo Set these to something more meaningfull!
3317  SELECT CASE(basis%INTERPOLATION_XI(1))
3318  CASE(basis_linear_simplex_interpolation)
3319  SELECT CASE(basis%NUMBER_OF_XI)
3320  CASE(1)
3321  basis%QUADRATURE%GAUSS_ORDER=2
3322  CASE(2)
3323  basis%QUADRATURE%GAUSS_ORDER=3
3324  CASE(3)
3325  basis%QUADRATURE%GAUSS_ORDER=3
3326  CASE DEFAULT
3327  local_error="The number of xi directions ("//trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))// &
3328  & ") is invalid"
3329  CALL flagerror(local_error,err,error,*999)
3330  END SELECT
3331  CASE(basis_quadratic_simplex_interpolation)
3332  SELECT CASE(basis%NUMBER_OF_XI)
3333  CASE(1)
3334  basis%QUADRATURE%GAUSS_ORDER=3
3335  CASE(2)
3336  basis%QUADRATURE%GAUSS_ORDER=3
3337  CASE(3)
3338  basis%QUADRATURE%GAUSS_ORDER=5
3339  CASE DEFAULT
3340  local_error="The number of xi directions ("//trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))// &
3341  & ") is invalid"
3342  CALL flagerror(local_error,err,error,*999)
3343  END SELECT
3344  CASE(basis_cubic_simplex_interpolation)
3345  SELECT CASE(basis%NUMBER_OF_XI)
3346  CASE(1)
3347  basis%QUADRATURE%GAUSS_ORDER=3
3348  CASE(2)
3349  basis%QUADRATURE%GAUSS_ORDER=5
3350  CASE(3)
3351  basis%QUADRATURE%GAUSS_ORDER=5
3352  CASE DEFAULT
3353  local_error="The number of xi directions ("//trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))// &
3354  & ") is invalid"
3355  CALL flagerror(local_error,err,error,*999)
3356  END SELECT
3357  CASE DEFAULT
3358  local_error="Interpolation xi value "//trim(number_to_vstring(basis%INTERPOLATION_XI(1),"*",err,error))// &
3359  & " in xi direction 1 is invalid"
3360  CALL flagerror(local_error,err,error,*999)
3361  END SELECT
3362  basis%QUADRATURE%NUMBER_OF_GAUSS_XI=basis%QUADRATURE%GAUSS_ORDER
3363  CASE DEFAULT
3364  local_error="Basis type value "//trim(number_to_vstring(basis%INTERPOLATION_XI(ni),"*",err,error))// &
3365  & " is invalid or not implemented"
3366  CALL flagerror(local_error,err,error,*999)
3367  END SELECT
3368  ENDIF
3369  ELSE
3370  CALL flagerror("Basis is not associated",err,error,*998)
3371  ENDIF
3372 
3373  exits("BASIS_QUADRATURE_INITIALISE")
3374  RETURN
3375 999 IF(ALLOCATED(basis%QUADRATURE%NUMBER_OF_GAUSS_XI)) DEALLOCATE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI)
3376 998 errorsexits("BASIS_QUADRATURE_INITIALISE",err,error)
3377  RETURN 1
3378 
3379  END SUBROUTINE basis_quadrature_initialise
3380 
3381  !
3382  !================================================================================================================================
3383  !
3384 
3386  SUBROUTINE basis_quadrature_number_of_gauss_xi_get(BASIS,QUADRATURE_NUMBER_OF_GAUSS_XI,ERR,ERROR,*)
3387 
3388  !Argument variables
3389  TYPE(basis_type), POINTER :: BASIS
3390  INTEGER(INTG), INTENT(OUT) :: QUADRATURE_NUMBER_OF_GAUSS_XI(:)
3391  INTEGER(INTG), INTENT(OUT) :: ERR
3392  TYPE(varying_string), INTENT(OUT) :: ERROR
3393  !Local Variables
3394  TYPE(varying_string) :: LOCAL_ERROR
3395 
3396  enters("BASIS_QUADRATURE_NUMBER_OF_GAUSS_XI_GET",err,error,*999)
3397 
3398  IF(ASSOCIATED(basis)) THEN
3399  IF(basis%BASIS_FINISHED) THEN
3400  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3401  IF(SIZE(quadrature_number_of_gauss_xi,1)>=SIZE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI,1)) THEN
3402  quadrature_number_of_gauss_xi=basis%QUADRATURE%NUMBER_OF_GAUSS_XI
3403  ELSE
3404  local_error="The size of QUADRATURE_NUMBER_OF_GAUSS_XI is too small. The supplied size is "// &
3405  & trim(number_to_vstring(SIZE(quadrature_number_of_gauss_xi,1),"*",err,error))//" and it needs to be >= "// &
3406  & trim(number_to_vstring(SIZE(basis%QUADRATURE%NUMBER_OF_GAUSS_XI,1),"*",err,error))//"."
3407  CALL flagerror(local_error,err,error,*999)
3408  ENDIF
3409  ELSE
3410  CALL flagerror("Quadrature basis is not associated.",err,error,*999)
3411  ENDIF
3412  ELSE
3413  CALL flagerror("Basis has not finished.",err,error,*999)
3414  ENDIF
3415  ELSE
3416  CALL flagerror("Basis is not associated.",err,error,*999)
3417  ENDIF
3418 
3419  exits("BASIS_QUADRATURE_NUMBER_OF_GAUSS_XI_GET")
3420  RETURN
3421 999 errorsexits("BASIS_QUADRATURE_NUMBER_OF_GAUSS_XI_GET",err,error)
3422  RETURN
3423  END SUBROUTINE basis_quadrature_number_of_gauss_xi_get
3424 
3425  !
3426  !================================================================================================================================
3427  !
3428  !================================================================================================================================
3429  !
3430 
3432  SUBROUTINE basis_quadrature_number_of_gauss_xi_set(BASIS,NUMBER_OF_GAUSS_XI,ERR,ERROR,*)
3433 
3434  !Argument variables
3435  TYPE(basis_type), POINTER :: BASIS
3436  INTEGER(INTG), INTENT(IN) :: NUMBER_OF_GAUSS_XI(:)
3437  INTEGER(INTG), INTENT(OUT) :: ERR
3438  TYPE(varying_string), INTENT(OUT) :: ERROR
3439  !Local Variables
3440  INTEGER(INTG) :: ni
3441  TYPE(varying_string) :: LOCAL_ERROR,LOCAL_WARNING
3442 
3443  enters("BASIS_QUADRATURE_NUMBER_OF_GAUSS_XI_SET",err,error,*999)
3444 
3445  IF(ASSOCIATED(basis)) THEN
3446  IF(basis%BASIS_FINISHED) THEN
3447  CALL flagerror("Basis has been finished.",err,error,*999)
3448  ELSE
3449  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3450  IF(SIZE(number_of_gauss_xi,1)==basis%NUMBER_OF_XI) THEN
3451  IF(any(number_of_gauss_xi<1)) CALL flagerror("Invalid number of gauss values.",err,error,*999)
3452  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI)=number_of_gauss_xi(1:basis%NUMBER_OF_XI)
3453  !Check the number of gauss points is sufficient for the interpolation order and flag a warning if not
3454  DO ni=1,basis%NUMBER_OF_XI
3455  SELECT CASE(basis%INTERPOLATION_XI(ni))
3456  CASE(basis_linear_lagrange_interpolation)
3457  IF(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)<2) THEN
3458  local_warning=trim(number_to_vstring(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),"*",err,error))// &
3459  & " Gauss points are insufficient for linear Lagrange interpolation."
3460  CALL flag_warning(local_warning,err,error,*999)
3461  ENDIF
3462  CASE(basis_quadratic_lagrange_interpolation)
3463  IF(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)<2) THEN
3464  local_warning=trim(number_to_vstring(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),"*",err,error))//&
3465  & " Gauss points are insufficient for quadratic Lagrange interpolation."
3466  CALL flag_warning(local_warning,err,error,*999)
3467  ENDIF
3468  CASE(basis_cubic_lagrange_interpolation)
3469  IF(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)<3) THEN
3470  local_warning=trim(number_to_vstring(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),"*",err,error))//&
3471  & " Gauss points are insufficient for cubic Lagrange interpolation."
3472  CALL flag_warning(local_warning,err,error,*999)
3473  ENDIF
3474  CASE(basis_quadratic1_hermite_interpolation,basis_quadratic2_hermite_interpolation)
3475  IF(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)<2) THEN
3476  local_warning=trim(number_to_vstring(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),"*",err,error))//&
3477  & " Gauss points are insufficient for quadratic Hermite interpolation."
3478  CALL flag_warning(local_warning,err,error,*999)
3479  ENDIF
3480  CASE(basis_cubic_hermite_interpolation)
3481  IF(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)<3) THEN
3482  local_warning=trim(number_to_vstring(basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni),"*",err,error))//&
3483  & " Gauss points are insufficient for cubic Hermite interpolation."
3484  CALL flag_warning(local_warning,err,error,*999)
3485  ENDIF
3486  CASE(basis_linear_simplex_interpolation)
3487  local_warning="For simplex elements please set quadrature order rather than number of gauss points."
3488  CALL flag_warning(local_warning,err,error,*999)
3489  CASE(basis_quadratic_simplex_interpolation)
3490  local_warning="For simplex elements please set quadrature order rather than number of gauss points."
3491  CALL flag_warning(local_warning,err,error,*999)
3492  CASE DEFAULT
3493  local_error="Interpolation xi value "//trim(number_to_vstring(basis%INTERPOLATION_XI(ni),"*",err,error))// &
3494  & " is invalid for xi direction "//trim(number_to_vstring(ni,"*",err,error))//"."
3495  CALL flagerror(local_error,err,error,*999)
3496  END SELECT
3497  ENDDO !xi
3498  ELSE
3499  local_error="The size of the number of Gauss array ("// &
3500  & trim(number_to_vstring(SIZE(number_of_gauss_xi,1),"*",err,error))// &
3501  & ") does not match the number of xi directions ("// &
3502  & trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))//") for basis number "// &
3503  & trim(number_to_vstring(basis%USER_NUMBER,"*",err,error))//"."
3504  CALL flagerror(local_error,err,error,*999)
3505  ENDIF
3506  ELSE
3507  CALL flagerror("Quadrature basis is not associated.",err,error,*999)
3508  ENDIF
3509  ENDIF
3510  ELSE
3511  CALL flagerror("Basis is not associated.",err,error,*999)
3512  ENDIF
3513 
3514  exits("BASIS_QUADRATURE_NUMBER_OF_GAUSS_XI_SET")
3515  RETURN
3516 999 errorsexits("BASIS_QUADRATURE_NUMBER_OF_GAUSS_XI_SET",err,error)
3517  RETURN 1
3518  END SUBROUTINE basis_quadrature_number_of_gauss_xi_set
3519  !
3520  !================================================================================================================================
3521  !
3522 
3524  SUBROUTINE basis_quadrature_single_gauss_xi_get(BASIS,SCHEME,GAUSS_POINT,GAUSS_XI,ERR,ERROR,*)
3525 
3526  !Argument variables
3527  TYPE(basis_type), POINTER :: BASIS
3528  INTEGER(INTG), INTENT(IN) :: SCHEME
3529  INTEGER(INTG), INTENT(IN) :: GAUSS_POINT
3530  REAL(DP), INTENT(OUT) :: GAUSS_XI(:)
3531  INTEGER(INTG), INTENT(OUT) :: ERR
3532  TYPE(varying_string), INTENT(OUT) :: ERROR
3533  !Local Variables
3534  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
3535  TYPE(varying_string) :: LOCAL_ERROR
3536 
3537  enters("BASIS_QUADRATURE_SINGLE_GAUSS_XI_GET",err,error,*999)
3538 
3539  IF(ASSOCIATED(basis)) THEN
3540  IF(basis%BASIS_FINISHED) THEN
3541  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3542  quadrature_scheme=>basis%QUADRATURE%QUADRATURE_SCHEME_MAP(scheme)%PTR
3543  IF(ASSOCIATED(quadrature_scheme)) THEN
3544  IF(SIZE(gauss_xi)==basis%NUMBER_OF_XI) THEN
3545  IF(gauss_point>0.AND.gauss_point<=quadrature_scheme%NUMBER_OF_GAUSS) THEN
3546  gauss_xi(:)=quadrature_scheme%GAUSS_POSITIONS(:,gauss_point)
3547  ELSE
3548  local_error="The specified Gauss point number of "// &
3549  & trim(number_to_vstring(gauss_point,"*",err,error))//" is invalid for the specified "// &
3550  & "quadrature scheme of the specified element for this field which has "// &
3551  & trim(number_to_vstring(quadrature_scheme%NUMBER_OF_GAUSS,"*",err,error))//" Gauss points."
3552  CALL flagerror(local_error,err,error,*999)
3553  ENDIF
3554  ELSE
3555  local_error="The number of xi values to return is invalid and needs to be "// &
3556  & trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))//" for the specified basis."
3557  CALL flagerror(local_error,err,error,*999)
3558  ENDIF
3559  ELSE
3560  CALL flagerror("The specified quadrature scheme is not associated for this basis.",err,error,*999)
3561  ENDIF
3562  ELSE
3563  CALL flagerror("Quadrature basis is not associated.",err,error,*999)
3564  ENDIF
3565  ELSE
3566  CALL flagerror("Basis has not finished.",err,error,*999)
3567  ENDIF
3568  ELSE
3569  CALL flagerror("Basis is not associated.",err,error,*999)
3570  ENDIF
3571 
3572  exits("BASIS_QUADRATURE_SINGLE_GAUSS_XI_GET")
3573  RETURN
3574 999 errorsexits("BASIS_QUADRATURE_SINGLE_GAUSS_XI_GET",err,error)
3575  RETURN
3576  END SUBROUTINE basis_quadrature_single_gauss_xi_get
3577 
3578  !
3579  !================================================================================================================================
3580  !
3581 
3583  SUBROUTINE basis_quadrature_multiple_gauss_xi_get(BASIS,SCHEME,GAUSS_POINTS,GAUSS_XI,ERR,ERROR,*)
3584 
3585  !Argument variables
3586  TYPE(basis_type), POINTER :: BASIS
3587  INTEGER(INTG), INTENT(IN) :: SCHEME
3588  INTEGER(INTG), INTENT(IN) :: GAUSS_POINTS(:)
3589  REAL(DP), INTENT(OUT) :: GAUSS_XI(:,:)
3590  INTEGER(INTG), INTENT(OUT) :: ERR
3591  TYPE(varying_string), INTENT(OUT) :: ERROR
3592  !Local Variables
3593  INTEGER(INTG) :: Gauss_point
3594  TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
3595  TYPE(varying_string) :: LOCAL_ERROR
3596 
3597  enters("BASIS_QUADRATURE_MULTIPLE_GAUSS_XI_GET",err,error,*999)
3598 
3599  IF(ASSOCIATED(basis)) THEN
3600  IF(basis%BASIS_FINISHED) THEN
3601  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3602  quadrature_scheme=>basis%QUADRATURE%QUADRATURE_SCHEME_MAP(scheme)%PTR
3603  IF(ASSOCIATED(quadrature_scheme)) THEN
3604  IF(SIZE(gauss_xi,1)==basis%NUMBER_OF_XI) THEN
3605  IF(SIZE(gauss_points)==0) THEN !Return all Gauss point xi locations.
3606  IF(SIZE(gauss_xi,2)==quadrature_scheme%NUMBER_OF_GAUSS) THEN
3607  gauss_xi=quadrature_scheme%GAUSS_POSITIONS
3608  ELSE
3609  local_error="The number of Gauss Points to return the xi values for is invalid and needs to be "// &
3610  & trim(number_to_vstring(quadrature_scheme%NUMBER_OF_GAUSS,"*",err,error))//"."
3611  CALL flagerror(local_error,err,error,*999)
3612  ENDIF
3613  ELSE !Return only specified Gauss point xi locations.
3614  DO gauss_point=1,SIZE(gauss_points)
3615  IF(gauss_points(gauss_point)>0.AND.gauss_points(gauss_point)<=quadrature_scheme%NUMBER_OF_GAUSS) THEN
3616  gauss_xi(:,gauss_point)=quadrature_scheme%GAUSS_POSITIONS(:,gauss_points(gauss_point))
3617  ELSE
3618  local_error="The specified Gauss point number of "// &
3619  & trim(number_to_vstring(gauss_points(gauss_point),"*",err,error))//" is invalid for the specified "// &
3620  & "quadrature scheme of the specified element for this field which has "// &
3621  & trim(number_to_vstring(quadrature_scheme%NUMBER_OF_GAUSS,"*",err,error))//" Gauss points."
3622  CALL flagerror(local_error,err,error,*999)
3623  ENDIF
3624  ENDDO
3625  ENDIF
3626  ELSE
3627  local_error="The number of xi values to return is invalid and needs to be "// &
3628  & trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))//" for the specified basis."
3629  CALL flagerror(local_error,err,error,*999)
3630  ENDIF
3631  ELSE
3632  CALL flagerror("The specified quadrature scheme is not associated for this basis.",err,error,*999)
3633  ENDIF
3634  ELSE
3635  CALL flagerror("Quadrature basis is not associated.",err,error,*999)
3636  ENDIF
3637  ELSE
3638  CALL flagerror("Basis has not finished.",err,error,*999)
3639  ENDIF
3640  ELSE
3641  CALL flagerror("Basis is not associated.",err,error,*999)
3642  ENDIF
3643 
3644  exits("BASIS_QUADRATURE_MULTIPLE_GAUSS_XI_GET")
3645  RETURN
3646 999 errorsexits("BASIS_QUADRATURE_MULTIPLE_GAUSS_XI_GET",err,error)
3647  RETURN
3648  END SUBROUTINE basis_quadrature_multiple_gauss_xi_get
3649 
3650  !
3651  !================================================================================================================================
3652  !
3653 
3655  SUBROUTINE basis_quadrature_order_get(BASIS,QUADRATURE_ORDER,ERR,ERROR,*)
3656 
3657  !Argument variables
3658  TYPE(basis_type), POINTER :: BASIS
3659  INTEGER(INTG), INTENT(OUT) :: QUADRATURE_ORDER
3660  INTEGER(INTG), INTENT(OUT) :: ERR
3661  TYPE(varying_string), INTENT(OUT) :: ERROR
3662  !Local Variables
3663 
3664  enters("BASIS_QUADRATURE_ORDER_GET",err,error,*999)
3665 
3666  IF(ASSOCIATED(basis)) THEN
3667  IF(basis%BASIS_FINISHED) THEN
3668  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3669  quadrature_order=basis%QUADRATURE%GAUSS_ORDER
3670  ELSE
3671  CALL flagerror("Quadrature basis is not associated.",err,error,*999)
3672  ENDIF
3673  ELSE
3674  CALL flagerror("Basis has not finished.",err,error,*999)
3675  ENDIF
3676  ELSE
3677  CALL flagerror("Basis is not associated.",err,error,*999)
3678  ENDIF
3679 
3680  exits("BASIS_QUADRATURE_ORDER_GET")
3681  RETURN
3682 999 errorsexits("BASIS_QUADRATURE_ORDER_GET",err,error)
3683  RETURN
3684  END SUBROUTINE basis_quadrature_order_get
3685 
3686  !
3687  !================================================================================================================================
3688  !
3689 
3691  SUBROUTINE basis_quadrature_order_set_number(USER_NUMBER,ORDER,ERR,ERROR,*)
3692 
3693  !Argument variables
3694  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
3695  INTEGER(INTG), INTENT(IN) :: ORDER
3696  INTEGER(INTG), INTENT(OUT) :: ERR
3697  TYPE(varying_string), INTENT(OUT) :: ERROR
3698  !Local Variables
3699  TYPE(basis_type), POINTER :: BASIS
3700 
3701  enters("BASIS_QUADRATURE_ORDER_SET_NUMBER",err,error,*999)
3702 
3703  CALL basis_user_number_find(user_number,basis,err,error,*999)
3704  CALL basis_quadrature_order_set(basis,order,err,error,*999)
3705 
3706  exits("BASIS_QUADRATURE_ORDER_SET_NUMBER")
3707  RETURN
3708 999 errorsexits("BASIS_QUADRATURE_ORDER_SET_NUMBER",err,error)
3709  RETURN 1
3710  END SUBROUTINE basis_quadrature_order_set_number
3711 
3712  !
3713  !================================================================================================================================
3714  !
3715 
3717  SUBROUTINE basis_quadrature_order_set_ptr(BASIS,ORDER,ERR,ERROR,*)
3718 
3719  !Argument variables
3720  TYPE(basis_type), POINTER :: BASIS
3721  INTEGER(INTG), INTENT(IN) :: ORDER
3722  INTEGER(INTG), INTENT(OUT) :: ERR
3723  TYPE(varying_string), INTENT(OUT) :: ERROR
3724  !Local Variables
3725  TYPE(varying_string) :: LOCAL_ERROR
3726 
3727  enters("BASIS_QUADRATURE_ORDER_SET_PTR",err,error,*999)
3728 
3729  IF(ASSOCIATED(basis)) THEN
3730  IF(basis%BASIS_FINISHED) THEN
3731  CALL flagerror("Basis has been finished",err,error,*999)
3732  ELSE
3733  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3734  IF(basis%TYPE==basis_simplex_type) THEN !Relax this i.e., use this to set gauss points in each direction for LHTP's???
3735  IF(order>1.AND.order<=5) THEN
3736  basis%QUADRATURE%GAUSS_ORDER=order
3737  ELSE
3738  local_error="An order value of "//trim(number_to_vstring(order,"*",err,error))// &
3739  & " is invalid. You must specify and order between 1 and 5."
3740  CALL flagerror(local_error,err,error,*999)
3741  ENDIF
3742  ELSE
3743  CALL flagerror("Can only set the quadrature order for simplex basis types.",err,error,*999)
3744  ENDIF
3745  ELSE
3746  CALL flagerror("Quadrature basis is not associated.",err,error,*999)
3747  ENDIF
3748  ENDIF
3749  ELSE
3750  CALL flagerror("Basis is not associated.",err,error,*999)
3751  ENDIF
3752 
3753  exits("BASIS_QUADRATURE_ORDER_SET_PTR")
3754  RETURN
3755 999 errorsexits("BASIS_QUADRATURE_ORDER_SET_PTR",err,error)
3756  RETURN 1
3757  END SUBROUTINE basis_quadrature_order_set_ptr
3758 
3759  !
3760  !================================================================================================================================
3761  !
3762 
3764  SUBROUTINE basis_quadrature_type_get(BASIS,QUADRATURE_TYPE,ERR,ERROR,*)
3765 
3766  !Argument variables
3767  TYPE(basis_type), POINTER :: BASIS
3768  INTEGER(INTG), INTENT(OUT) :: QUADRATURE_TYPE
3769  INTEGER(INTG), INTENT(OUT) :: ERR
3770  TYPE(varying_string), INTENT(OUT) :: ERROR
3771  !Local Variables
3772 
3773  enters("BASIS_QUADRATURE_TYPE_GET",err,error,*999)
3774 
3775  IF(ASSOCIATED(basis)) THEN
3776  IF(basis%BASIS_FINISHED) THEN
3777  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3778  quadrature_type=basis%QUADRATURE%TYPE
3779  ELSE
3780  CALL flagerror("Basis quadrature basis is not associated.",err,error,*999)
3781  ENDIF
3782  ELSE
3783  CALL flagerror("Basis has not finished.",err,error,*999)
3784  ENDIF
3785  ELSE
3786  CALL flagerror("Basis is not associated.",err,error,*999)
3787  ENDIF
3788 
3789  exits("BASIS_QUADRATURE_TYPE_GET")
3790  RETURN
3791 999 errorsexits("BASIS_QUADRATURE_TYPE_GET",err,error)
3792  RETURN
3793  END SUBROUTINE basis_quadrature_type_get
3794 
3795  !
3796  !================================================================================================================================
3797  !
3798 
3800  SUBROUTINE basis_quadrature_type_set_number(USER_NUMBER,TYPE,ERR,ERROR,*)
3801 
3802  !Argument variables
3803  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
3804  INTEGER(INTG), INTENT(IN) :: TYPE
3805  INTEGER(INTG), INTENT(OUT) :: ERR
3806  TYPE(varying_string), INTENT(OUT) :: ERROR
3807  !Local Variables
3808  TYPE(basis_type), POINTER :: BASIS
3809 
3810  enters("BASIS_QUADRATURE_TYPE_SET_NUMBER",err,error,*999)
3811 
3812  CALL basis_user_number_find(user_number,basis,err,error,*999)
3813  CALL basis_quadrature_type_set_ptr(basis,TYPE,ERR,ERROR,*999)
3814 
3815  exits("BASIS_QUADRATURE_TYPE_SET_NUMBER")
3816  RETURN
3817 999 errorsexits("BASIS_QUADRATURE_TYPE_SET_NUMBER",err,error)
3818  RETURN 1
3819  END SUBROUTINE basis_quadrature_type_set_number
3820 
3821  !
3822  !================================================================================================================================
3823  !
3824 
3826  SUBROUTINE basis_quadrature_type_set_ptr(BASIS,TYPE,ERR,ERROR,*)
3827 
3828  !Argument variables
3829  TYPE(basis_type), POINTER :: BASIS
3830  INTEGER(INTG), INTENT(IN) :: TYPE
3831  INTEGER(INTG), INTENT(OUT) :: ERR
3832  TYPE(varying_string), INTENT(OUT) :: ERROR
3833  !Local Variables
3834  TYPE(varying_string) :: LOCAL_ERROR
3835 
3836  enters("BASIS_QUADRATURE_TYPE_SET_PTR",err,error,*999)
3837 
3838  IF(ASSOCIATED(basis)) THEN
3839  IF(basis%BASIS_FINISHED) THEN
3840  CALL flagerror("Basis has been finished.",err,error,*999)
3841  ELSE
3842  IF(ASSOCIATED(basis%QUADRATURE%BASIS)) THEN
3843  SELECT CASE(type)
3844  CASE(basis_gauss_legendre_quadrature)
3845  basis%QUADRATURE%TYPE=basis_gauss_legendre_quadrature
3846  CASE(basis_gauss_laguerre_quadrature)
3847  basis%QUADRATURE%TYPE=basis_gauss_laguerre_quadrature
3848  CALL flagerror("Gauss Laguerre quadrature is not implemented.",err,error,*999)
3849  CASE(basis_guass_hermite_quadrature)
3850  basis%QUADRATURE%TYPE=basis_guass_hermite_quadrature
3851  CALL flagerror("Gauss Hermite quadrature is not implemented.",err,error,*999)
3852  CASE DEFAULT
3853  local_error="Quadrature type "//trim(number_to_vstring(TYPE,"*",ERR,ERROR))//" is invalid."
3854  CALL flagerror(local_error,err,error,*999)
3855  END SELECT
3856  ELSE
3857  CALL flagerror("Basis quadrature basis is not associated.",err,error,*999)
3858  ENDIF
3859  ENDIF
3860  ELSE
3861  CALL flagerror("Basis is not associated.",err,error,*999)
3862  ENDIF
3863 
3864  exits("BASIS_QUADRATURE_TYPE_SET_PTR")
3865  RETURN
3866 999 errorsexits("BASIS_QUADRATURE_TYPE_SET_PTR",err,error)
3867  RETURN 1
3868  END SUBROUTINE basis_quadrature_type_set_ptr
3869 
3870  !
3871  !================================================================================================================================
3872  !
3873 
3875  SUBROUTINE basis_quadraturelocalfacegaussevaluateset(BASIS,FACE_GAUSS_EVALUATE,ERR,ERROR,*)
3876 
3877  !Argument variables
3878  TYPE(basis_type), POINTER :: BASIS
3879  LOGICAL, INTENT(IN) :: FACE_GAUSS_EVALUATE
3880  INTEGER(INTG), INTENT(OUT) :: ERR
3881  TYPE(varying_string), INTENT(OUT) :: ERROR
3882 
3883  enters("Basis_QuadratureLocalFaceGaussEvaluateSet",err,error,*999)
3884 
3885  IF(ASSOCIATED(basis)) THEN
3886  IF(basis%BASIS_FINISHED) THEN
3887  CALL flagerror("Basis has been finished.",err,error,*999)
3888  ELSE
3889  basis%QUADRATURE%EVALUATE_FACE_GAUSS=face_gauss_evaluate
3890  ENDIF
3891  ELSE
3892  CALL flagerror("Basis is not associated.",err,error,*999)
3893  ENDIF
3894 
3895  exits("Basis_QuadratureLocalFaceGaussEvaluateSet")
3896  RETURN
3897 999 errorsexits("Basis_QuadratureLocalFaceGaussEvaluateSet",err,error)
3898  RETURN 1
3899 
3900  END SUBROUTINE basis_quadraturelocalfacegaussevaluateset
3901 
3902  !
3903  !================================================================================================================================
3904  !
3905 
3908  SUBROUTINE basis_simplex_basis_create(BASIS,ERR,ERROR,*)
3909 
3910  !Argument variables
3911  TYPE(basis_type), POINTER :: BASIS
3912  INTEGER(INTG), INTENT(OUT) :: ERR
3913  TYPE(varying_string), INTENT(OUT) :: ERROR
3914  !Local Variables
3915  INTEGER(INTG) :: MAX_NUM_NODES,ni,nn,ns
3916  INTEGER(INTG), ALLOCATABLE :: NODES_IN_FACE(:)
3917  TYPE(varying_string) :: LOCAL_ERROR
3918 
3919  enters("BASIS_SIMPLEX_BASIS_CREATE",err,error,*999)
3920 
3921  IF(ASSOCIATED(basis)) THEN
3922  IF(basis%TYPE==basis_simplex_type) THEN
3923  basis%NUMBER_OF_XI_COORDINATES=basis%NUMBER_OF_XI+1 !Simplex bases have an additional area coordinate
3924  ALLOCATE(basis%INTERPOLATION_TYPE(basis%NUMBER_OF_XI_COORDINATES),stat=err)
3925  IF(err/=0) CALL flagerror("Could not allocate INTERPOLATION_TYPE array.",err,error,*999)
3926  ALLOCATE(basis%INTERPOLATION_ORDER(basis%NUMBER_OF_XI_COORDINATES),stat=err)
3927  IF(err/=0) CALL flagerror("Could not allocate INTERPOLATION_ORDER array.",err,error,*999)
3928  ALLOCATE(basis%NUMBER_OF_NODES_XIC(basis%NUMBER_OF_XI_COORDINATES),stat=err)
3929  IF(err/=0) CALL flagerror("Could not allocate NUMBER_OF_NODES_XIC array.",err,error,*999)
3930  basis%DEGENERATE=.false.
3931  basis%NUMBER_OF_COLLAPSED_XI=0
3932  SELECT CASE(basis%NUMBER_OF_XI)
3933  CASE(1)
3934  basis%NUMBER_OF_PARTIAL_DERIVATIVES=3
3935  SELECT CASE(basis%INTERPOLATION_XI(1))
3936  CASE(basis_linear_simplex_interpolation)
3937  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
3938  basis%INTERPOLATION_ORDER(1)=basis_linear_interpolation_order
3939  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
3940  basis%INTERPOLATION_ORDER(2)=basis_linear_interpolation_order
3941  basis%NUMBER_OF_NODES_XIC(1)=2
3942  basis%NUMBER_OF_NODES_XIC(2)=2
3943  max_num_nodes=2
3944  basis%NUMBER_OF_NODES=2
3945  CASE(basis_quadratic_simplex_interpolation)
3946  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
3947  basis%INTERPOLATION_ORDER(1)=basis_quadratic_interpolation_order
3948  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
3949  basis%INTERPOLATION_ORDER(2)=basis_quadratic_interpolation_order
3950  basis%NUMBER_OF_NODES_XIC(1)=3
3951  basis%NUMBER_OF_NODES_XIC(2)=3
3952  max_num_nodes=3
3953  basis%NUMBER_OF_NODES=3
3954  CASE(basis_cubic_simplex_interpolation)
3955  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
3956  basis%INTERPOLATION_ORDER(1)=basis_cubic_interpolation_order
3957  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
3958  basis%INTERPOLATION_ORDER(2)=basis_cubic_interpolation_order
3959  basis%NUMBER_OF_NODES_XIC(1)=4
3960  basis%NUMBER_OF_NODES_XIC(2)=4
3961  max_num_nodes=4
3962  basis%NUMBER_OF_NODES=4
3963  CASE DEFAULT
3964  CALL flagerror("Invalid interpolation type.",err,error,*999)
3965  END SELECT
3966  CASE(2)
3967  basis%NUMBER_OF_PARTIAL_DERIVATIVES=6
3968  SELECT CASE(basis%INTERPOLATION_XI(2))
3969  CASE(basis_linear_simplex_interpolation)
3970  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
3971  basis%INTERPOLATION_ORDER(1)=basis_linear_interpolation_order
3972  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
3973  basis%INTERPOLATION_ORDER(2)=basis_linear_interpolation_order
3974  basis%INTERPOLATION_TYPE(3)=basis_simplex_interpolation
3975  basis%INTERPOLATION_ORDER(3)=basis_linear_interpolation_order
3976  basis%NUMBER_OF_NODES_XIC(1)=2
3977  basis%NUMBER_OF_NODES_XIC(2)=2
3978  basis%NUMBER_OF_NODES_XIC(3)=2
3979  max_num_nodes=2
3980  basis%NUMBER_OF_NODES=3
3981  CASE(basis_quadratic_simplex_interpolation)
3982  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
3983  basis%INTERPOLATION_ORDER(1)=basis_quadratic_interpolation_order
3984  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
3985  basis%INTERPOLATION_ORDER(2)=basis_quadratic_interpolation_order
3986  basis%INTERPOLATION_TYPE(3)=basis_simplex_interpolation
3987  basis%INTERPOLATION_ORDER(3)=basis_quadratic_interpolation_order
3988  basis%NUMBER_OF_NODES_XIC(1)=3
3989  basis%NUMBER_OF_NODES_XIC(2)=3
3990  basis%NUMBER_OF_NODES_XIC(3)=3
3991  max_num_nodes=3
3992  basis%NUMBER_OF_NODES=6
3993  CASE(basis_cubic_simplex_interpolation)
3994  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
3995  basis%INTERPOLATION_ORDER(1)=basis_cubic_interpolation_order
3996  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
3997  basis%INTERPOLATION_ORDER(2)=basis_cubic_interpolation_order
3998  basis%INTERPOLATION_TYPE(3)=basis_simplex_interpolation
3999  basis%INTERPOLATION_ORDER(3)=basis_cubic_interpolation_order
4000  basis%NUMBER_OF_NODES_XIC(1)=4
4001  basis%NUMBER_OF_NODES_XIC(2)=4
4002  basis%NUMBER_OF_NODES_XIC(3)=4
4003  max_num_nodes=4
4004  basis%NUMBER_OF_NODES=10
4005  CASE DEFAULT
4006  CALL flagerror("Invalid interpolation type.",err,error,*999)
4007  END SELECT
4008  CASE(3)
4009  basis%NUMBER_OF_PARTIAL_DERIVATIVES=11
4010  SELECT CASE(basis%INTERPOLATION_XI(3))
4011  CASE(basis_linear_simplex_interpolation)
4012  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
4013  basis%INTERPOLATION_ORDER(1)=basis_linear_interpolation_order
4014  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
4015  basis%INTERPOLATION_ORDER(2)=basis_linear_interpolation_order
4016  basis%INTERPOLATION_TYPE(3)=basis_simplex_interpolation
4017  basis%INTERPOLATION_ORDER(3)=basis_linear_interpolation_order
4018  basis%INTERPOLATION_TYPE(4)=basis_simplex_interpolation
4019  basis%INTERPOLATION_ORDER(4)=basis_linear_interpolation_order
4020  basis%NUMBER_OF_NODES_XIC(1)=2
4021  basis%NUMBER_OF_NODES_XIC(2)=2
4022  basis%NUMBER_OF_NODES_XIC(3)=2
4023  basis%NUMBER_OF_NODES_XIC(4)=2
4024  max_num_nodes=2
4025  basis%NUMBER_OF_NODES=4
4026  CASE(basis_quadratic_simplex_interpolation)
4027  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
4028  basis%INTERPOLATION_ORDER(1)=basis_quadratic_interpolation_order
4029  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
4030  basis%INTERPOLATION_ORDER(2)=basis_quadratic_interpolation_order
4031  basis%INTERPOLATION_TYPE(3)=basis_simplex_interpolation
4032  basis%INTERPOLATION_ORDER(3)=basis_quadratic_interpolation_order
4033  basis%INTERPOLATION_TYPE(4)=basis_simplex_interpolation
4034  basis%INTERPOLATION_ORDER(4)=basis_quadratic_interpolation_order
4035  basis%NUMBER_OF_NODES_XIC(1)=3
4036  basis%NUMBER_OF_NODES_XIC(2)=3
4037  basis%NUMBER_OF_NODES_XIC(3)=3
4038  basis%NUMBER_OF_NODES_XIC(4)=3
4039  max_num_nodes=3
4040  basis%NUMBER_OF_NODES=10
4041  CASE(basis_cubic_simplex_interpolation)
4042  basis%INTERPOLATION_TYPE(1)=basis_simplex_interpolation
4043  basis%INTERPOLATION_ORDER(1)=basis_cubic_interpolation_order
4044  basis%INTERPOLATION_TYPE(2)=basis_simplex_interpolation
4045  basis%INTERPOLATION_ORDER(2)=basis_cubic_interpolation_order
4046  basis%INTERPOLATION_TYPE(3)=basis_simplex_interpolation
4047  basis%INTERPOLATION_ORDER(3)=basis_cubic_interpolation_order
4048  basis%INTERPOLATION_TYPE(4)=basis_simplex_interpolation
4049  basis%INTERPOLATION_ORDER(4)=basis_cubic_interpolation_order
4050  basis%NUMBER_OF_NODES_XIC(1)=4
4051  basis%NUMBER_OF_NODES_XIC(2)=4
4052  basis%NUMBER_OF_NODES_XIC(3)=4
4053  basis%NUMBER_OF_NODES_XIC(4)=4
4054  max_num_nodes=4
4055  basis%NUMBER_OF_NODES=20
4056  CASE DEFAULT
4057  CALL flagerror("Invalid interpolation type.",err,error,*999)
4058  END SELECT
4059  CASE DEFAULT
4060  CALL flagerror("Invalid number of xi directions.",err,error,*999)
4061  END SELECT
4062 
4063  ALLOCATE(basis%NODE_AT_COLLAPSE(basis%NUMBER_OF_NODES),stat=err)
4064  IF(err/=0) CALL flagerror("Could not allocate node at collapse.",err,error,*999)
4065  basis%NODE_AT_COLLAPSE=.false.
4066 
4067  ALLOCATE(basis%NODE_POSITION_INDEX(basis%NUMBER_OF_NODES,basis%NUMBER_OF_XI_COORDINATES),stat=err)
4068  IF(err/=0) CALL flagerror("Could not allocate NODE_POSITION_INDEX.",err,error,*999)
4069  SELECT CASE(basis%NUMBER_OF_XI_COORDINATES)
4070  CASE(2)
4071  ALLOCATE(basis%NODE_POSITION_INDEX_INV(max_num_nodes,max_num_nodes,1,1),stat=err)
4072  CASE(3)
4073  ALLOCATE(basis%NODE_POSITION_INDEX_INV(max_num_nodes,max_num_nodes,max_num_nodes,1),stat=err)
4074  CASE(4)
4075  ALLOCATE(basis%NODE_POSITION_INDEX_INV(max_num_nodes,max_num_nodes,max_num_nodes,max_num_nodes),stat=err)
4076  CASE DEFAULT
4077  CALL flagerror("Invalid number of coordinates.",err,error,*999)
4078  END SELECT
4079  IF(err/=0) CALL flagerror("Could not allocate NODE_POSITION_INDEX_INV.",err,error,*999)
4080  basis%NODE_POSITION_INDEX_INV=0
4081 
4082  !Determine the node position index and it's inverse
4083  SELECT CASE(basis%NUMBER_OF_XI)
4084  CASE(1)
4085  SELECT CASE(basis%INTERPOLATION_ORDER(1))
4086  CASE(basis_linear_interpolation_order)
4087  !Node 1
4088  basis%NODE_POSITION_INDEX(1,1)=2
4089  basis%NODE_POSITION_INDEX(1,2)=1
4090  basis%NODE_POSITION_INDEX_INV(2,1,1,1)=1
4091  !Node 2
4092  basis%NODE_POSITION_INDEX(2,1)=1
4093  basis%NODE_POSITION_INDEX(2,2)=2
4094  basis%NODE_POSITION_INDEX_INV(1,2,1,1)=2
4095  CASE(basis_quadratic_interpolation_order)
4096  !Node 1
4097  basis%NODE_POSITION_INDEX(1,1)=3
4098  basis%NODE_POSITION_INDEX(1,2)=1
4099  basis%NODE_POSITION_INDEX_INV(3,1,1,1)=1
4100  !Node 2
4101  basis%NODE_POSITION_INDEX(2,1)=2
4102  basis%NODE_POSITION_INDEX(2,2)=2
4103  basis%NODE_POSITION_INDEX_INV(2,2,1,1)=2
4104  !Node 3
4105  basis%NODE_POSITION_INDEX(3,1)=1
4106  basis%NODE_POSITION_INDEX(3,2)=3
4107  basis%NODE_POSITION_INDEX_INV(1,3,1,1)=3
4108  CASE(basis_cubic_interpolation_order)
4109  !Node 1
4110  basis%NODE_POSITION_INDEX(1,1)=4
4111  basis%NODE_POSITION_INDEX(1,2)=1
4112  basis%NODE_POSITION_INDEX_INV(4,1,1,1)=1
4113  !Node 2
4114  basis%NODE_POSITION_INDEX(2,1)=3
4115  basis%NODE_POSITION_INDEX(2,2)=2
4116  basis%NODE_POSITION_INDEX_INV(3,2,1,1)=2
4117  !Node 3
4118  basis%NODE_POSITION_INDEX(3,1)=2
4119  basis%NODE_POSITION_INDEX(3,2)=3
4120  basis%NODE_POSITION_INDEX_INV(2,3,1,1)=3
4121  !Node 4
4122  basis%NODE_POSITION_INDEX(4,1)=1
4123  basis%NODE_POSITION_INDEX(4,2)=4
4124  basis%NODE_POSITION_INDEX_INV(1,4,1,1)=4
4125  CASE DEFAULT
4126  CALL flagerror("Invalid interpolation order.",err,error,*999)
4127  END SELECT
4128  CASE(2)
4129  SELECT CASE(basis%INTERPOLATION_ORDER(1))
4130  CASE(basis_linear_interpolation_order)
4131  !Node 1
4132  basis%NODE_POSITION_INDEX(1,1)=2
4133  basis%NODE_POSITION_INDEX(1,2)=1
4134  basis%NODE_POSITION_INDEX(1,3)=1
4135  basis%NODE_POSITION_INDEX_INV(2,1,1,1)=1
4136  !Node 2
4137  basis%NODE_POSITION_INDEX(2,1)=1
4138  basis%NODE_POSITION_INDEX(2,2)=2
4139  basis%NODE_POSITION_INDEX(2,3)=1
4140  basis%NODE_POSITION_INDEX_INV(1,2,1,1)=2
4141  !Node 3
4142  basis%NODE_POSITION_INDEX(3,1)=1
4143  basis%NODE_POSITION_INDEX(3,2)=1
4144  basis%NODE_POSITION_INDEX(3,3)=2
4145  basis%NODE_POSITION_INDEX_INV(1,1,2,1)=3
4146  CASE(basis_quadratic_interpolation_order)
4147  !Node 1
4148  basis%NODE_POSITION_INDEX(1,1)=3
4149  basis%NODE_POSITION_INDEX(1,2)=1
4150  basis%NODE_POSITION_INDEX(1,3)=1
4151  basis%NODE_POSITION_INDEX_INV(3,1,1,1)=1
4152  !Node 2
4153  basis%NODE_POSITION_INDEX(2,1)=1
4154  basis%NODE_POSITION_INDEX(2,2)=3
4155  basis%NODE_POSITION_INDEX(2,3)=1
4156  basis%NODE_POSITION_INDEX_INV(1,3,1,1)=2
4157  !Node 3
4158  basis%NODE_POSITION_INDEX(3,1)=1
4159  basis%NODE_POSITION_INDEX(3,2)=1
4160  basis%NODE_POSITION_INDEX(3,3)=3
4161  basis%NODE_POSITION_INDEX_INV(1,1,3,1)=3
4162  !Node 4
4163  basis%NODE_POSITION_INDEX(4,1)=2
4164  basis%NODE_POSITION_INDEX(4,2)=2
4165  basis%NODE_POSITION_INDEX(4,3)=1
4166  basis%NODE_POSITION_INDEX_INV(2,2,1,1)=4
4167  !Node 5
4168  basis%NODE_POSITION_INDEX(5,1)=1
4169  basis%NODE_POSITION_INDEX(5,2)=2
4170  basis%NODE_POSITION_INDEX(5,3)=2
4171  basis%NODE_POSITION_INDEX_INV(1,2,2,1)=5
4172  !Node 6
4173  basis%NODE_POSITION_INDEX(6,1)=2
4174  basis%NODE_POSITION_INDEX(6,2)=1
4175  basis%NODE_POSITION_INDEX(6,3)=2
4176  basis%NODE_POSITION_INDEX_INV(2,1,2,1)=6
4177  CASE(basis_cubic_interpolation_order)
4178  !Node 1
4179  basis%NODE_POSITION_INDEX(1,1)=4
4180  basis%NODE_POSITION_INDEX(1,2)=1
4181  basis%NODE_POSITION_INDEX(1,3)=1
4182  basis%NODE_POSITION_INDEX_INV(4,1,1,1)=1
4183  !Node 2
4184  basis%NODE_POSITION_INDEX(2,1)=1
4185  basis%NODE_POSITION_INDEX(2,2)=4
4186  basis%NODE_POSITION_INDEX(2,3)=1
4187  basis%NODE_POSITION_INDEX_INV(1,4,1,1)=2
4188  !Node 3
4189  basis%NODE_POSITION_INDEX(3,1)=1
4190  basis%NODE_POSITION_INDEX(3,2)=1
4191  basis%NODE_POSITION_INDEX(3,3)=4
4192  basis%NODE_POSITION_INDEX_INV(1,1,4,1)=3
4193  !Node 4
4194  basis%NODE_POSITION_INDEX(4,1)=3
4195  basis%NODE_POSITION_INDEX(4,2)=2
4196  basis%NODE_POSITION_INDEX(4,3)=1
4197  basis%NODE_POSITION_INDEX_INV(3,2,1,1)=4
4198  !Node 5
4199  basis%NODE_POSITION_INDEX(5,1)=2
4200  basis%NODE_POSITION_INDEX(5,2)=3
4201  basis%NODE_POSITION_INDEX(5,3)=1
4202  basis%NODE_POSITION_INDEX_INV(2,3,1,1)=5
4203  !Node 6
4204  basis%NODE_POSITION_INDEX(6,1)=1
4205  basis%NODE_POSITION_INDEX(6,2)=3
4206  basis%NODE_POSITION_INDEX(6,3)=2
4207  basis%NODE_POSITION_INDEX_INV(1,3,2,1)=6
4208  !Node 7
4209  basis%NODE_POSITION_INDEX(7,1)=1
4210  basis%NODE_POSITION_INDEX(7,2)=2
4211  basis%NODE_POSITION_INDEX(7,3)=3
4212  basis%NODE_POSITION_INDEX_INV(1,2,3,1)=7
4213  !Node 8
4214  basis%NODE_POSITION_INDEX(8,1)=2
4215  basis%NODE_POSITION_INDEX(8,2)=1
4216  basis%NODE_POSITION_INDEX(8,3)=3
4217  basis%NODE_POSITION_INDEX_INV(2,1,3,1)=8
4218  !Node 9
4219  basis%NODE_POSITION_INDEX(9,1)=3
4220  basis%NODE_POSITION_INDEX(9,2)=1
4221  basis%NODE_POSITION_INDEX(9,3)=2
4222  basis%NODE_POSITION_INDEX_INV(3,1,2,1)=9
4223  !Node 10
4224  basis%NODE_POSITION_INDEX(10,1)=2
4225  basis%NODE_POSITION_INDEX(10,2)=2
4226  basis%NODE_POSITION_INDEX(10,3)=2
4227  basis%NODE_POSITION_INDEX_INV(2,2,2,1)=10
4228  CASE DEFAULT
4229  CALL flagerror("Invalid interpolation order",err,error,*999)
4230  END SELECT
4231  CASE(3)
4232  SELECT CASE(basis%INTERPOLATION_ORDER(1))
4233  CASE(basis_linear_interpolation_order)
4234  !Node 1
4235  basis%NODE_POSITION_INDEX(1,1)=2
4236  basis%NODE_POSITION_INDEX(1,2)=1
4237  basis%NODE_POSITION_INDEX(1,3)=1
4238  basis%NODE_POSITION_INDEX(1,4)=1
4239  basis%NODE_POSITION_INDEX_INV(2,1,1,1)=1
4240  !Node 2
4241  basis%NODE_POSITION_INDEX(2,1)=1
4242  basis%NODE_POSITION_INDEX(2,2)=2
4243  basis%NODE_POSITION_INDEX(2,3)=1
4244  basis%NODE_POSITION_INDEX(2,4)=1
4245  basis%NODE_POSITION_INDEX_INV(1,2,1,1)=2
4246  !Node 3
4247  basis%NODE_POSITION_INDEX(3,1)=1
4248  basis%NODE_POSITION_INDEX(3,2)=1
4249  basis%NODE_POSITION_INDEX(3,3)=2
4250  basis%NODE_POSITION_INDEX(3,4)=1
4251  basis%NODE_POSITION_INDEX_INV(1,1,2,1)=3
4252  !Node 4
4253  basis%NODE_POSITION_INDEX(4,1)=1
4254  basis%NODE_POSITION_INDEX(4,2)=1
4255  basis%NODE_POSITION_INDEX(4,3)=1
4256  basis%NODE_POSITION_INDEX(4,4)=2
4257  basis%NODE_POSITION_INDEX_INV(1,1,1,2)=4
4258 
4259  ALLOCATE(nodes_in_face(12),stat=err)
4260  IF(err/=0) CALL flagerror("Could not allocate NODES_IN_FACE.",err,error,*999)
4261  nodes_in_face(:)=[2,3,4,1,3,4,1,2,4,1,2,3] !12 Nodes
4262 
4263  CASE(basis_quadratic_interpolation_order)
4264  !Node 1
4265  basis%NODE_POSITION_INDEX(1,1)=3
4266  basis%NODE_POSITION_INDEX(1,2)=1
4267  basis%NODE_POSITION_INDEX(1,3)=1
4268  basis%NODE_POSITION_INDEX(1,4)=1
4269  basis%NODE_POSITION_INDEX_INV(3,1,1,1)=1
4270  !Node 2
4271  basis%NODE_POSITION_INDEX(2,1)=1
4272  basis%NODE_POSITION_INDEX(2,2)=3
4273  basis%NODE_POSITION_INDEX(2,3)=1
4274  basis%NODE_POSITION_INDEX(2,4)=1
4275  basis%NODE_POSITION_INDEX_INV(1,3,1,1)=2
4276  !Node 3
4277  basis%NODE_POSITION_INDEX(3,1)=1
4278  basis%NODE_POSITION_INDEX(3,2)=1
4279  basis%NODE_POSITION_INDEX(3,3)=3
4280  basis%NODE_POSITION_INDEX(3,4)=1
4281  basis%NODE_POSITION_INDEX_INV(1,1,3,1)=3
4282  !Node 4
4283  basis%NODE_POSITION_INDEX(4,1)=1
4284  basis%NODE_POSITION_INDEX(4,2)=1
4285  basis%NODE_POSITION_INDEX(4,3)=1
4286  basis%NODE_POSITION_INDEX(4,4)=3
4287  basis%NODE_POSITION_INDEX_INV(1,1,1,3)=4
4288  !Node 5
4289  basis%NODE_POSITION_INDEX(5,1)=2
4290  basis%NODE_POSITION_INDEX(5,2)=2
4291  basis%NODE_POSITION_INDEX(5,3)=1
4292  basis%NODE_POSITION_INDEX(5,4)=1
4293  basis%NODE_POSITION_INDEX_INV(2,2,1,1)=5
4294  !Node 6
4295  basis%NODE_POSITION_INDEX(6,1)=2
4296  basis%NODE_POSITION_INDEX(6,2)=1
4297  basis%NODE_POSITION_INDEX(6,3)=2
4298  basis%NODE_POSITION_INDEX(6,4)=1
4299  basis%NODE_POSITION_INDEX_INV(2,1,2,1)=6
4300  !Node 7
4301  basis%NODE_POSITION_INDEX(7,1)=2
4302  basis%NODE_POSITION_INDEX(7,2)=1
4303  basis%NODE_POSITION_INDEX(7,3)=1
4304  basis%NODE_POSITION_INDEX(7,4)=2
4305  basis%NODE_POSITION_INDEX_INV(2,1,1,2)=7
4306  !Node 8
4307  basis%NODE_POSITION_INDEX(8,1)=1
4308  basis%NODE_POSITION_INDEX(8,2)=2
4309  basis%NODE_POSITION_INDEX(8,3)=2
4310  basis%NODE_POSITION_INDEX(8,4)=1
4311  basis%NODE_POSITION_INDEX_INV(1,2,2,1)=8
4312  !Node 9
4313  basis%NODE_POSITION_INDEX(9,1)=1
4314  basis%NODE_POSITION_INDEX(9,2)=1
4315  basis%NODE_POSITION_INDEX(9,3)=2
4316  basis%NODE_POSITION_INDEX(9,4)=2
4317  basis%NODE_POSITION_INDEX_INV(1,1,2,2)=9
4318  !Node 10
4319  basis%NODE_POSITION_INDEX(10,1)=1
4320  basis%NODE_POSITION_INDEX(10,2)=2
4321  basis%NODE_POSITION_INDEX(10,3)=1
4322  basis%NODE_POSITION_INDEX(10,4)=2
4323  basis%NODE_POSITION_INDEX_INV(1,2,1,2)=10
4324 
4325  ALLOCATE(nodes_in_face(24),stat=err)
4326  IF(err/=0) CALL flagerror("Could not allocate NODES_IN_FACE.",err,error,*999)
4327  nodes_in_face(:)=[2,3,4,8,9,10,1,3,4,6,9,7,1,2,4,5,10,7,1,2,3,5,8,6] !24 Nodes
4328 
4329  CASE(basis_cubic_interpolation_order)
4330  !Node 1
4331  basis%NODE_POSITION_INDEX(1,1)=4
4332  basis%NODE_POSITION_INDEX(1,2)=1
4333  basis%NODE_POSITION_INDEX(1,3)=1
4334  basis%NODE_POSITION_INDEX(1,4)=1
4335  basis%NODE_POSITION_INDEX_INV(4,1,1,1)=1
4336  !Node 2
4337  basis%NODE_POSITION_INDEX(2,1)=1
4338  basis%NODE_POSITION_INDEX(2,2)=4
4339  basis%NODE_POSITION_INDEX(2,3)=1
4340  basis%NODE_POSITION_INDEX(2,4)=1
4341  basis%NODE_POSITION_INDEX_INV(1,4,1,1)=2
4342  !Node 3
4343  basis%NODE_POSITION_INDEX(3,1)=1
4344  basis%NODE_POSITION_INDEX(3,2)=1
4345  basis%NODE_POSITION_INDEX(3,3)=4
4346  basis%NODE_POSITION_INDEX(3,4)=1
4347  basis%NODE_POSITION_INDEX_INV(1,1,4,1)=3
4348  !Node 4
4349  basis%NODE_POSITION_INDEX(4,1)=1
4350  basis%NODE_POSITION_INDEX(4,2)=1
4351  basis%NODE_POSITION_INDEX(4,3)=1
4352  basis%NODE_POSITION_INDEX(4,4)=4
4353  basis%NODE_POSITION_INDEX_INV(1,1,1,4)=4
4354  !Node 5
4355  basis%NODE_POSITION_INDEX(5,1)=3
4356  basis%NODE_POSITION_INDEX(5,2)=2
4357  basis%NODE_POSITION_INDEX(5,3)=1
4358  basis%NODE_POSITION_INDEX(5,4)=1
4359  basis%NODE_POSITION_INDEX_INV(3,2,1,1)=5
4360  !Node 6
4361  basis%NODE_POSITION_INDEX(6,1)=2
4362  basis%NODE_POSITION_INDEX(6,2)=3
4363  basis%NODE_POSITION_INDEX(6,3)=1
4364  basis%NODE_POSITION_INDEX(6,4)=1
4365  basis%NODE_POSITION_INDEX_INV(2,3,1,1)=6
4366  !Node 7
4367  basis%NODE_POSITION_INDEX(7,1)=3
4368  basis%NODE_POSITION_INDEX(7,2)=1
4369  basis%NODE_POSITION_INDEX(7,3)=2
4370  basis%NODE_POSITION_INDEX(7,4)=1
4371  basis%NODE_POSITION_INDEX_INV(3,1,2,1)=7
4372  !Node 8
4373  basis%NODE_POSITION_INDEX(8,1)=2
4374  basis%NODE_POSITION_INDEX(8,2)=1
4375  basis%NODE_POSITION_INDEX(8,3)=3
4376  basis%NODE_POSITION_INDEX(8,4)=1
4377  basis%NODE_POSITION_INDEX_INV(2,1,3,1)=8
4378  !Node 9
4379  basis%NODE_POSITION_INDEX(9,1)=3
4380  basis%NODE_POSITION_INDEX(9,2)=1
4381  basis%NODE_POSITION_INDEX(9,3)=1
4382  basis%NODE_POSITION_INDEX(9,4)=2
4383  basis%NODE_POSITION_INDEX_INV(3,1,1,2)=9
4384  !Node 10
4385  basis%NODE_POSITION_INDEX(10,1)=2
4386  basis%NODE_POSITION_INDEX(10,2)=1
4387  basis%NODE_POSITION_INDEX(10,3)=1
4388  basis%NODE_POSITION_INDEX(10,4)=3
4389  basis%NODE_POSITION_INDEX_INV(2,1,1,3)=10
4390  !Node 11
4391  basis%NODE_POSITION_INDEX(11,1)=1
4392  basis%NODE_POSITION_INDEX(11,2)=3
4393  basis%NODE_POSITION_INDEX(11,3)=2
4394  basis%NODE_POSITION_INDEX(11,4)=1
4395  basis%NODE_POSITION_INDEX_INV(1,3,2,1)=11
4396  !Node 12
4397  basis%NODE_POSITION_INDEX(12,1)=1
4398  basis%NODE_POSITION_INDEX(12,2)=2
4399  basis%NODE_POSITION_INDEX(12,3)=3
4400  basis%NODE_POSITION_INDEX(12,4)=1
4401  basis%NODE_POSITION_INDEX_INV(1,2,3,1)=12
4402  !Node 13
4403  basis%NODE_POSITION_INDEX(13,1)=1
4404  basis%NODE_POSITION_INDEX(13,2)=1
4405  basis%NODE_POSITION_INDEX(13,3)=3
4406  basis%NODE_POSITION_INDEX(13,4)=2
4407  basis%NODE_POSITION_INDEX_INV(1,1,3,2)=13
4408  !Node 14
4409  basis%NODE_POSITION_INDEX(14,1)=1
4410  basis%NODE_POSITION_INDEX(14,2)=1
4411  basis%NODE_POSITION_INDEX(14,3)=2
4412  basis%NODE_POSITION_INDEX(14,4)=3
4413  basis%NODE_POSITION_INDEX_INV(1,1,2,3)=14
4414  !Node 15
4415  basis%NODE_POSITION_INDEX(15,1)=1
4416  basis%NODE_POSITION_INDEX(15,2)=3
4417  basis%NODE_POSITION_INDEX(15,3)=1
4418  basis%NODE_POSITION_INDEX(15,4)=2
4419  basis%NODE_POSITION_INDEX_INV(1,3,1,2)=15
4420  !Node 16
4421  basis%NODE_POSITION_INDEX(16,1)=1
4422  basis%NODE_POSITION_INDEX(16,2)=2
4423  basis%NODE_POSITION_INDEX(16,3)=1
4424  basis%NODE_POSITION_INDEX(16,4)=3
4425  basis%NODE_POSITION_INDEX_INV(1,2,1,3)=16
4426  !Node 17
4427  basis%NODE_POSITION_INDEX(17,1)=2
4428  basis%NODE_POSITION_INDEX(17,2)=2
4429  basis%NODE_POSITION_INDEX(17,3)=2
4430  basis%NODE_POSITION_INDEX(17,4)=1
4431  basis%NODE_POSITION_INDEX_INV(2,2,2,1)=17
4432  !Node 18
4433  basis%NODE_POSITION_INDEX(18,1)=2
4434  basis%NODE_POSITION_INDEX(18,2)=2
4435  basis%NODE_POSITION_INDEX(18,3)=1
4436  basis%NODE_POSITION_INDEX(18,4)=2
4437  basis%NODE_POSITION_INDEX_INV(2,2,1,2)=18
4438  !Node 19
4439  basis%NODE_POSITION_INDEX(19,1)=2
4440  basis%NODE_POSITION_INDEX(19,2)=1
4441  basis%NODE_POSITION_INDEX(19,3)=2
4442  basis%NODE_POSITION_INDEX(19,4)=2
4443  basis%NODE_POSITION_INDEX_INV(2,1,2,2)=19
4444  !Node 20
4445  basis%NODE_POSITION_INDEX(20,1)=1
4446  basis%NODE_POSITION_INDEX(20,2)=2
4447  basis%NODE_POSITION_INDEX(20,3)=2
4448  basis%NODE_POSITION_INDEX(20,4)=2
4449  basis%NODE_POSITION_INDEX_INV(1,2,2,2)=20
4450 
4451  ALLOCATE(nodes_in_face(40),stat=err)
4452  IF(err/=0) CALL flagerror("Could not allocate NODES_IN_FACE.",err,error,*999)
4453  nodes_in_face(:)=[2,3,4,11,12,13,14,16,15,20,1,3,4,7,8,13,14,10,9,&
4454  &19,1,2,4,5,6,15,16,10,9,18,1,2,3,5,6,14,12,8,7,17] !40 nodes
4455 
4456  CASE DEFAULT
4457  CALL flagerror("Invalid interpolation order.",err,error,*999)
4458  END SELECT
4459  CASE DEFAULT
4460  CALL flagerror("Invalid number of xi directions.",err,error,*999)
4461  END SELECT
4462  !Calculate the maximum number of derivatives (1 for simplex bases) and the number of element parameters
4463  basis%MAXIMUM_NUMBER_OF_DERIVATIVES=1
4464  basis%NUMBER_OF_ELEMENT_PARAMETERS=basis%NUMBER_OF_NODES
4465  !Now set up the number of derivatives and derivative order index
4466  ALLOCATE(basis%NUMBER_OF_DERIVATIVES(basis%NUMBER_OF_NODES),stat=err)
4467  IF(err/=0) CALL flagerror("Could not allocate NUMBER_OF_DERIVATIVES.",err,error,*999)
4468  ALLOCATE(basis%DERIVATIVE_ORDER_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES,basis%NUMBER_OF_XI), &
4469  & stat=err)
4470  IF(err/=0) CALL flagerror("Could not allocate DERIVATIVE_ORDER_INDEX.",err,error,*999)
4471  ALLOCATE(basis%DERIVATIVE_ORDER_INDEX_INV(first_part_deriv,first_part_deriv,first_part_deriv,basis%NUMBER_OF_NODES), &
4472  & stat=err)
4473  IF(err/=0) CALL flagerror("Could not allocate DERIVATIVE_ORDER_INDEX_INV.",err,error,*999)
4474  ALLOCATE(basis%PARTIAL_DERIVATIVE_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES),stat=err)
4475  IF(err/=0) CALL flagerror("Could not allocate PARTIAL_DERIVATIVE_INDEX.",err,error,*999)
4476  ALLOCATE(basis%ELEMENT_PARAMETER_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES),stat=err)
4477  IF(err/=0) CALL flagerror("Could not allocate ELEMENT_PARAMETER_INDEX.",err,error,*999)
4478  ALLOCATE(basis%ELEMENT_PARAMETER_INDEX_INV(2,basis%NUMBER_OF_ELEMENT_PARAMETERS),stat=err)
4479  IF(err/=0) CALL flagerror("Could not allocate ELEMENT_PARAMETER_INDEX_INV.",err,error,*999)
4480  !Set the derivative order index and its inverse, the element parameter index and the partial derivative index.
4481  ns=0
4482  basis%DERIVATIVE_ORDER_INDEX_INV=0
4483  DO nn=1,basis%NUMBER_OF_NODES
4484  basis%NUMBER_OF_DERIVATIVES(nn)=1
4485  DO ni=1,basis%NUMBER_OF_XI
4486  basis%DERIVATIVE_ORDER_INDEX(1,nn,ni)=1
4487  ENDDO !ni
4488  ns=ns+1
4489  basis%ELEMENT_PARAMETER_INDEX(1,nn)=ns
4490  basis%ELEMENT_PARAMETER_INDEX_INV(1,ns)=nn
4491  basis%ELEMENT_PARAMETER_INDEX_INV(2,ns)=1
4492  basis%PARTIAL_DERIVATIVE_INDEX(1,nn)=no_part_deriv
4493  basis%DERIVATIVE_ORDER_INDEX_INV(basis%DERIVATIVE_ORDER_INDEX(1,nn,1),1,1,nn)=1
4494  ENDDO !nn
4495 
4496  !Set up the line and face information
4497  SELECT CASE(basis%NUMBER_OF_XI)
4498  CASE(1)
4499  basis%NUMBER_OF_LOCAL_LINES=1
4500  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4501  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local line.",err,error,*999)
4502  ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4503  IF(err/=0) CALL flagerror("Could not allocate local line xi direction.",err,error,*999)
4504  basis%LOCAL_LINE_XI_DIRECTION(1)=1
4505  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1),basis%NUMBER_OF_LOCAL_LINES),stat=err)
4506  IF(err/=0) CALL flagerror("Could not allocate node numbers in local line.",err,error,*999)
4507  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1),basis%NUMBER_OF_LOCAL_LINES),stat=err)
4508  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local line.",err,error,*999)
4509  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=no_part_deriv
4510  ALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1)**2,basis%NUMBER_OF_LOCAL_LINES) &
4511  & ,stat=err)
4512  IF(err/=0) CALL flagerror("Could not allocate element parameters in local line.",err,error,*999)
4513  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE=1
4514  !Set the line values
4515  SELECT CASE(basis%INTERPOLATION_ORDER(1))
4516  CASE(basis_linear_interpolation_order)
4517  !Line 1
4518  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=2
4519  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4520  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=2
4521  CASE(basis_quadratic_interpolation_order)
4522  !Line 1
4523  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=3
4524  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4525  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=2
4526  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,1)=3
4527  CASE(basis_cubic_interpolation_order)
4528  !Line 1
4529  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=4
4530  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4531  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=2
4532  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,1)=3
4533  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,1)=4
4534  CASE DEFAULT
4535  local_error="Interpolation order "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(1),"*", &
4536  & err,error))//" is invalid for a simplex basis type."
4537  CALL flagerror(local_error,err,error,*999)
4538  END SELECT
4539  CASE(2)
4540  !Allocate and calculate the lines
4541  !Simplex hence three local lines
4542  basis%NUMBER_OF_LOCAL_LINES=3
4543  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4544  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local line.",err,error,*999)
4545  ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4546  IF(err/=0) CALL flagerror("Could not allocate local line xi direction",err,error,*999)
4547  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),basis%NUMBER_OF_LOCAL_LINES),stat=err)
4548  IF(err/=0) CALL flagerror("Could not allocate node numbers in local line.",err,error,*999)
4549  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),basis%NUMBER_OF_LOCAL_LINES),stat=err)
4550  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local line.",err,error,*999)
4551  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=no_part_deriv
4552  ALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC)**2,basis%NUMBER_OF_LOCAL_LINES) &
4553  & ,stat=err)
4554  IF(err/=0) CALL flagerror("Could not allocate element parameters in local line.",err,error,*999)
4555  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE=1
4556  ALLOCATE(basis%LOCAL_XI_NORMAL(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4557  IF(err/=0) CALL flagerror("Could not allocate local line normal.",err,error,*999)
4558  !Set the line values
4559  SELECT CASE(basis%INTERPOLATION_ORDER(1))
4560  CASE(basis_linear_interpolation_order)
4561  !Line 1
4562  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=2
4563  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4564  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=2
4565  basis%LOCAL_LINE_XI_DIRECTION(1)=1
4566  basis%LOCAL_XI_NORMAL(1)=3
4567  !Line 2
4568  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(2)=2
4569  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,2)=1
4570  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,2)=3
4571  basis%LOCAL_LINE_XI_DIRECTION(2)=2
4572  basis%LOCAL_XI_NORMAL(2)=2
4573  !Line 3
4574  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(3)=2
4575  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,3)=2
4576  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,3)=3
4577  basis%LOCAL_LINE_XI_DIRECTION(3)=3
4578  basis%LOCAL_XI_NORMAL(3)=1
4579  CASE(basis_quadratic_interpolation_order)
4580  !Line 1
4581  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=3
4582  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4583  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=4
4584  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,1)=2
4585  basis%LOCAL_LINE_XI_DIRECTION(1)=1
4586  basis%LOCAL_XI_NORMAL(1)=3
4587  !Line 2
4588  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(2)=3
4589  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,2)=1
4590  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,2)=6
4591  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,2)=3
4592  basis%LOCAL_LINE_XI_DIRECTION(2)=2
4593  basis%LOCAL_XI_NORMAL(2)=2
4594  !Line 3
4595  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(3)=3
4596  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,3)=2
4597  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,3)=5
4598  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,3)=3
4599  basis%LOCAL_LINE_XI_DIRECTION(3)=3
4600  basis%LOCAL_XI_NORMAL(3)=1
4601  CASE(basis_cubic_interpolation_order)
4602  !Line 1
4603  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=4
4604  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4605  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=4
4606  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,1)=5
4607  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,1)=2
4608  basis%LOCAL_LINE_XI_DIRECTION(1)=1
4609  basis%LOCAL_XI_NORMAL(1)=3
4610  !Line 2
4611  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(2)=4
4612  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,2)=1
4613  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,2)=9
4614  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,2)=8
4615  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,2)=3
4616  basis%LOCAL_LINE_XI_DIRECTION(2)=2
4617  basis%LOCAL_XI_NORMAL(2)=2
4618  !Line 3
4619  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(3)=4
4620  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,3)=2
4621  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,3)=6
4622  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,3)=7
4623  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,3)=3
4624  basis%LOCAL_LINE_XI_DIRECTION(3)=3
4625  basis%LOCAL_XI_NORMAL(3)=1
4626  CASE DEFAULT
4627  local_error="Interpolation order "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(1),"*",err,error))// &
4628  & " is invalid for a simplex basis type."
4629  CALL flagerror(local_error,err,error,*999)
4630  END SELECT
4631  CASE(3)
4632  basis%NUMBER_OF_LOCAL_LINES=6
4633  basis%NUMBER_OF_LOCAL_FACES=4
4634 
4635  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4636  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local line.",err,error,*999)
4637  ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_FACE(basis%NUMBER_OF_LOCAL_FACES),stat=err)
4638  IF(err/=0) CALL flagerror("Could not allocate number of nodes in local face.",err,error,*999)
4639 
4640  ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4641  IF(err/=0) CALL flagerror("Could not allocate local line xi direction.",err,error,*999)
4642  ALLOCATE(basis%LOCAL_FACE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_FACES),stat=err)
4643  IF(err/=0) CALL flagerror("Could not allocate local face xi direction.",err,error,*999)
4644 
4645  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),basis%NUMBER_OF_LOCAL_LINES),stat=err)
4646  IF(err/=0) CALL flagerror("Could not allocate node numbers in local line.",err,error,*999)
4647  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC),basis%NUMBER_OF_LOCAL_LINES),stat=err)
4648  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local line.",err,error,*999)
4649  basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=no_part_deriv
4650  ALLOCATE(basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE(maxval(basis%NUMBER_OF_NODES_XIC)**2,basis%NUMBER_OF_LOCAL_LINES), &
4651  & stat=err)
4652  IF(err/=0) CALL flagerror("Could not allocate element parameters in local line.",err,error,*999)
4653  basis%ELEMENT_PARAMETERS_IN_LOCAL_LINE=1
4654 
4655  SELECT CASE(basis%INTERPOLATION_ORDER(1))
4656  CASE(basis_linear_interpolation_order)
4657  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_FACE(3,basis%NUMBER_OF_LOCAL_FACES),stat=err)
4658  IF(err/=0) CALL flagerror("Could not allocate node numbers in local face.",err,error,*999)
4659  !\todo Number of local face node derivatives currenlty set to 1 (calculation of BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE for simplex elements has not been implemented yet)
4660  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:1,3,basis%NUMBER_OF_LOCAL_FACES),stat=err)
4661  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local face.",err,error,*999)
4662  CASE(basis_quadratic_interpolation_order)
4663  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_FACE(6,basis%NUMBER_OF_LOCAL_FACES),stat=err)
4664  IF(err/=0) CALL flagerror("Could not allocate node numbers in local face.",err,error,*999)
4665  !\todo Number of local face node derivatives currenlty set to 1 (calculation of BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE for simplex elements has not been implemented yet)
4666  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:1,6,basis%NUMBER_OF_LOCAL_FACES),stat=err)
4667  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local face.",err,error,*999)
4668  CASE(basis_cubic_interpolation_order)
4669  ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_FACE(10,basis%NUMBER_OF_LOCAL_FACES),stat=err)
4670  IF(err/=0) CALL flagerror("Could not allocate node numbers in local face.",err,error,*999)
4671  !\todo Number of local face node derivatives currenlty set to 1 (calculation of BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE for simplex elements has not been implemented yet)
4672  ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:1,10,basis%NUMBER_OF_LOCAL_FACES),stat=err)
4673  IF(err/=0) CALL flagerror("Could not allocate derivative numbers in local face.",err,error,*999)
4674  CASE DEFAULT
4675  local_error="Interpolation order "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(1),"*",err,error))// &
4676  & " is invalid for a simplex basis type."
4677  CALL flagerror(local_error,err,error,*999)
4678  END SELECT
4679  basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(1,:,:)=no_part_deriv
4680  basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0,:,:)=1
4681 
4682  ALLOCATE(basis%LOCAL_XI_NORMAL(basis%NUMBER_OF_LOCAL_LINES),stat=err)
4683  IF(err/=0) CALL flagerror("Could not allocate local line normal.",err,error,*999)
4684 
4685  !Set the line values
4686  SELECT CASE(basis%INTERPOLATION_ORDER(1))
4687  CASE(basis_linear_interpolation_order)
4688  !Line 1
4689  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=2
4690  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4691  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=2
4692  basis%LOCAL_LINE_XI_DIRECTION(1)=1
4693  !Line 2
4694  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(2)=2
4695  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,2)=1
4696  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,2)=3
4697  basis%LOCAL_LINE_XI_DIRECTION(2)=1
4698  !Line 3
4699  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(3)=2
4700  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,3)=1
4701  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,3)=4
4702  basis%LOCAL_LINE_XI_DIRECTION(3)=1
4703  !Line 4
4704  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(4)=2
4705  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,4)=2
4706  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,4)=3
4707  basis%LOCAL_LINE_XI_DIRECTION(4)=2
4708  !Line 5
4709  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(5)=2
4710  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,5)=2
4711  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,5)=4
4712  basis%LOCAL_LINE_XI_DIRECTION(5)=2
4713  !Line 6
4714  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(6)=2
4715  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,6)=3
4716  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,6)=4
4717  basis%LOCAL_LINE_XI_DIRECTION(6)=3
4718  !Face 1
4719  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(1)=3
4720  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,1)=2
4721  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,1)=3
4722  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,1)=4
4723  basis%LOCAL_FACE_XI_DIRECTION(1)=1
4724  basis%LOCAL_XI_NORMAL(1)=1
4725  !Face 2
4726  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(2)=3
4727  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,2)=1
4728  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,2)=4
4729  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,2)=3
4730  basis%LOCAL_FACE_XI_DIRECTION(2)=2
4731  basis%LOCAL_XI_NORMAL(2)=2
4732  !Face 3
4733  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(3)=3
4734  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,3)=1
4735  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,3)=2
4736  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,3)=4
4737  basis%LOCAL_FACE_XI_DIRECTION(3)=3
4738  basis%LOCAL_XI_NORMAL(3)=3
4739  !Face 4
4740  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(4)=3
4741  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,4)=1
4742  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,4)=3
4743  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,4)=2
4744  basis%LOCAL_FACE_XI_DIRECTION(4)=4
4745  basis%LOCAL_XI_NORMAL(4)=4
4746  CASE(basis_quadratic_interpolation_order)
4747  !Line 1
4748  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=3
4749  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4750  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=5
4751  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,1)=2
4752  basis%LOCAL_LINE_XI_DIRECTION(1)=1
4753  !Line 2
4754  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(2)=3
4755  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,2)=1
4756  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,2)=6
4757  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,2)=3
4758  basis%LOCAL_LINE_XI_DIRECTION(2)=1
4759  !Line 3
4760  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(3)=3
4761  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,3)=1
4762  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,3)=7
4763  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,3)=4
4764  basis%LOCAL_LINE_XI_DIRECTION(3)=1
4765  !Line 4
4766  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(4)=3
4767  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,4)=2
4768  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,4)=8
4769  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,4)=3
4770  basis%LOCAL_LINE_XI_DIRECTION(4)=2
4771  !Line 5
4772  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(5)=3
4773  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,5)=2
4774  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,5)=10
4775  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,5)=4
4776  basis%LOCAL_LINE_XI_DIRECTION(5)=2
4777  !Line 6
4778  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(6)=3
4779  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,6)=3
4780  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,6)=9
4781  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,6)=4
4782  basis%LOCAL_LINE_XI_DIRECTION(6)=3
4783  !Face 1
4784  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(1)=6
4785  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,1)=2
4786  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,1)=3
4787  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,1)=4
4788  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,1)=8
4789  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,1)=9
4790  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,1)=10
4791  basis%LOCAL_FACE_XI_DIRECTION(1)=1
4792  basis%LOCAL_XI_NORMAL(1)=1
4793  !Face 2
4794  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(2)=6
4795  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,2)=1
4796  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,2)=4
4797  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,2)=3
4798  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,2)=7
4799  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,2)=9
4800  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,2)=6
4801  basis%LOCAL_FACE_XI_DIRECTION(2)=2
4802  basis%LOCAL_XI_NORMAL(2)=2
4803  !Face 3
4804  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(3)=6
4805  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,3)=1
4806  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,3)=2
4807  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,3)=4
4808  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,3)=5
4809  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,3)=10
4810  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,3)=7
4811  basis%LOCAL_FACE_XI_DIRECTION(3)=3
4812  basis%LOCAL_XI_NORMAL(3)=3
4813  !Face 4
4814  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(4)=6
4815  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,4)=1
4816  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,4)=3
4817  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,4)=2
4818  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,4)=6
4819  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,4)=8
4820  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,4)=5
4821  basis%LOCAL_FACE_XI_DIRECTION(4)=4
4822  basis%LOCAL_XI_NORMAL(4)=4
4823  CASE(basis_cubic_interpolation_order)
4824  !Line 1
4825  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=4
4826  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,1)=1
4827  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,1)=5
4828  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,1)=6
4829  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,1)=2
4830  basis%LOCAL_LINE_XI_DIRECTION(1)=1
4831  !Line 2
4832  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(2)=4
4833  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,2)=1
4834  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,2)=7
4835  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,2)=8
4836  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,2)=3
4837  basis%LOCAL_LINE_XI_DIRECTION(2)=1
4838  !Line 3
4839  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(3)=4
4840  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,3)=1
4841  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,3)=9
4842  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,3)=10
4843  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,3)=4
4844  basis%LOCAL_LINE_XI_DIRECTION(3)=1
4845  !Line 4
4846  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(4)=4
4847  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,4)=2
4848  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,4)=11
4849  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,4)=12
4850  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,4)=3
4851  basis%LOCAL_LINE_XI_DIRECTION(4)=2
4852  !Line 5
4853  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(5)=4
4854  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,5)=2
4855  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,5)=15
4856  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,5)=16
4857  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,5)=4
4858  basis%LOCAL_LINE_XI_DIRECTION(5)=2
4859  !Line 6
4860  basis%NUMBER_OF_NODES_IN_LOCAL_LINE(6)=4
4861  basis%NODE_NUMBERS_IN_LOCAL_LINE(1,6)=3
4862  basis%NODE_NUMBERS_IN_LOCAL_LINE(2,6)=13
4863  basis%NODE_NUMBERS_IN_LOCAL_LINE(3,6)=14
4864  basis%NODE_NUMBERS_IN_LOCAL_LINE(4,6)=4
4865  basis%LOCAL_LINE_XI_DIRECTION(6)=3
4866  !Face 1
4867  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(1)=10
4868  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,1)=2
4869  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,1)=3
4870  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,1)=4
4871  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,1)=11
4872  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,1)=12
4873  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,1)=13
4874  basis%NODE_NUMBERS_IN_LOCAL_FACE(7,1)=14
4875  basis%NODE_NUMBERS_IN_LOCAL_FACE(8,1)=16
4876  basis%NODE_NUMBERS_IN_LOCAL_FACE(9,1)=15
4877  basis%NODE_NUMBERS_IN_LOCAL_FACE(10,1)=20
4878  basis%LOCAL_FACE_XI_DIRECTION(1)=1
4879  basis%LOCAL_XI_NORMAL(1)=1
4880  !Face 2
4881  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(2)=10
4882  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,2)=1
4883  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,2)=4
4884  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,2)=3
4885  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,2)=9
4886  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,2)=10
4887  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,2)=14
4888  basis%NODE_NUMBERS_IN_LOCAL_FACE(7,2)=13
4889  basis%NODE_NUMBERS_IN_LOCAL_FACE(8,2)=8
4890  basis%NODE_NUMBERS_IN_LOCAL_FACE(9,2)=7
4891  basis%NODE_NUMBERS_IN_LOCAL_FACE(10,2)=19
4892  basis%LOCAL_FACE_XI_DIRECTION(2)=2
4893  basis%LOCAL_XI_NORMAL(2)=2
4894  !Face 3
4895  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(3)=10
4896  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,3)=1
4897  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,3)=2
4898  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,3)=4
4899  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,3)=5
4900  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,3)=6
4901  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,3)=15
4902  basis%NODE_NUMBERS_IN_LOCAL_FACE(7,3)=16
4903  basis%NODE_NUMBERS_IN_LOCAL_FACE(8,3)=10
4904  basis%NODE_NUMBERS_IN_LOCAL_FACE(9,3)=9
4905  basis%NODE_NUMBERS_IN_LOCAL_FACE(10,3)=18
4906  basis%LOCAL_FACE_XI_DIRECTION(3)=3
4907  basis%LOCAL_XI_NORMAL(3)=3
4908  !Face 4
4909  basis%NUMBER_OF_NODES_IN_LOCAL_FACE(4)=10
4910  basis%NODE_NUMBERS_IN_LOCAL_FACE(1,4)=1
4911  basis%NODE_NUMBERS_IN_LOCAL_FACE(2,4)=3
4912  basis%NODE_NUMBERS_IN_LOCAL_FACE(3,4)=2
4913  basis%NODE_NUMBERS_IN_LOCAL_FACE(4,4)=7
4914  basis%NODE_NUMBERS_IN_LOCAL_FACE(5,4)=8
4915  basis%NODE_NUMBERS_IN_LOCAL_FACE(6,4)=12
4916  basis%NODE_NUMBERS_IN_LOCAL_FACE(7,4)=11
4917  basis%NODE_NUMBERS_IN_LOCAL_FACE(8,4)=6
4918  basis%NODE_NUMBERS_IN_LOCAL_FACE(9,4)=5
4919  basis%NODE_NUMBERS_IN_LOCAL_FACE(10,4)=17
4920  basis%LOCAL_FACE_XI_DIRECTION(4)=4
4921  basis%LOCAL_XI_NORMAL(4)=4
4922  CASE DEFAULT
4923  local_error="Interpolation order "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(1),"*",err,error))// &
4924  & " is invalid for a simplex basis type."
4925  CALL flagerror(local_error,err,error,*999)
4926  END SELECT
4927  CASE DEFAULT
4928  CALL flagerror("Invalid number of xi directions.",err,error,*999)
4929  END SELECT
4930 
4931  CALL basis_quadrature_create(basis,err,error,*999)
4932 
4933  ELSE
4934  CALL flagerror("Basis is not a simplex basis.",err,error,*999)
4935  ENDIF
4936  ELSE
4937  CALL flagerror("Basis is not associated.",err,error,*999)
4938  ENDIF
4939 
4940  exits("BASIS_SIMPLEX_BASIS_CREATE")
4941  RETURN
4942 999 errorsexits("BASIS_SIMPLEX_BASIS_CREATE",err,error)
4943  RETURN 1
4944  END SUBROUTINE basis_simplex_basis_create
4945 
4946  !
4947  !================================================================================================================================
4948  !
4949 
4952  SUBROUTINE basis_simplex_family_create(BASIS,ERR,ERROR,*)
4953 
4954  !Argument variables
4955  TYPE(basis_type), POINTER :: BASIS
4956  INTEGER(INTG), INTENT(OUT) :: ERR
4957  TYPE(varying_string), INTENT(OUT) :: ERROR
4958  !Local Variables
4959  INTEGER(INTG) :: DUMMY_ERR,ni,ni2,FACE_XI(2),FACE_XI2(2)
4960  LOGICAL :: LINE_BASIS_DONE,FACE_BASIS_DONE
4961  TYPE(basis_type), POINTER :: NEW_SUB_BASIS
4962  TYPE(varying_string) :: DUMMY_ERROR
4963 
4964  NULLIFY(new_sub_basis)
4965 
4966  enters("BASIS_SIMPLEX_FAMILY_CREATE",err,error,*999)
4967 
4968  IF(ASSOCIATED(basis)) THEN
4969  !Create the main (parent) basis
4970  CALL basis_simplex_basis_create(basis,err,error,*999)
4971  IF(basis%NUMBER_OF_XI>1) THEN
4972  !Create the line bases as sub-basis types
4973 ! ALLOCATE(BASIS%LINE_BASES(BASIS%NUMBER_OF_XI),STAT=ERR)
4974  IF (basis%NUMBER_OF_XI .EQ. 2) THEN
4975  ALLOCATE(basis%LINE_BASES(basis%NUMBER_OF_XI+1),stat=err)
4976  ELSE
4977  ALLOCATE(basis%LINE_BASES(basis%NUMBER_OF_XI),stat=err)
4978  ENDIF
4979 
4980  IF(err/=0) CALL flagerror("Could not allocate basis line bases.",err,error,*999)
4981  DO ni=1,basis%NUMBER_OF_XI
4982  line_basis_done=.false.
4983  NULLIFY(new_sub_basis)
4984  DO ni2=1,ni-1
4985  IF(basis%INTERPOLATION_XI(ni2)==basis%INTERPOLATION_XI(ni).AND. &
4986  basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni2)==basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)) THEN
4987  line_basis_done=.true.
4988  EXIT
4989  ENDIF
4990  ENDDO !ni2
4991  IF(line_basis_done) THEN
4992  basis%LINE_BASES(ni)%PTR=>basis%LINE_BASES(ni2)%PTR
4993  ELSE
4994  !Create the new sub-basis
4995  CALL basis_sub_basis_create(basis,1,[ni],new_sub_basis,err,error,*999)
4996  !Fill in the basis information
4997  CALL basis_simplex_basis_create(new_sub_basis,err,error,*999)
4998  basis%LINE_BASES(ni)%PTR=>new_sub_basis
4999  ENDIF
5000  ENDDO !ni
5001 
5002  IF (basis%NUMBER_OF_XI .EQ. 2) THEN
5003  basis%LINE_BASES(basis%NUMBER_OF_XI+1)%PTR=>basis%LINE_BASES(basis%NUMBER_OF_XI)%PTR
5004  ENDIF
5005 
5006  IF(basis%NUMBER_OF_XI>2) THEN
5007  !Set up face basis functions
5008  ALLOCATE(basis%FACE_BASES(basis%NUMBER_OF_XI),stat=err)
5009  IF(err/=0) CALL flagerror("Could not allocate basis face bases.",err,error,*999)
5010  DO ni=1,basis%NUMBER_OF_XI
5011  !Determine the face xi directions that lie in this xi direction
5012  face_xi(1)=other_xi_directions3(ni,2,1)
5013  face_xi(2)=other_xi_directions3(ni,3,1)
5014  face_basis_done=.false.
5015  NULLIFY(new_sub_basis)
5016  DO ni2=1,ni-1
5017  face_xi2(1)=other_xi_directions3(ni2,2,1)
5018  face_xi2(2)=other_xi_directions3(ni2,3,1)
5019  IF(basis%INTERPOLATION_XI(face_xi2(1))==basis%INTERPOLATION_XI(face_xi(1)).AND. &
5020  & basis%INTERPOLATION_XI(face_xi2(2))==basis%INTERPOLATION_XI(face_xi(2)).AND. &
5021  & basis%QUADRATURE%NUMBER_OF_GAUSS_XI(face_xi2(1))==basis%QUADRATURE%NUMBER_OF_GAUSS_XI(face_xi(1)).AND. &
5022  & basis%QUADRATURE%NUMBER_OF_GAUSS_XI(face_xi2(2))==basis%QUADRATURE%NUMBER_OF_GAUSS_XI(face_xi(1))) THEN
5023  face_basis_done=.true.
5024  EXIT
5025  ENDIF
5026  ENDDO !ni2
5027  IF(face_basis_done) THEN
5028  basis%FACE_BASES(ni)%PTR=>basis%FACE_BASES(ni2)%PTR
5029  ELSE
5030  !Create the new sub-basis
5031  CALL basis_sub_basis_create(basis,2,[face_xi(1),face_xi(2)],new_sub_basis,err,error,*999)
5032  !Fill in the basis information
5033  CALL basis_simplex_basis_create(new_sub_basis,err,error,*999)
5034  new_sub_basis%LINE_BASES(1)%PTR=>basis%LINE_BASES(face_xi(1))%PTR
5035  new_sub_basis%LINE_BASES(2)%PTR=>basis%LINE_BASES(face_xi(2))%PTR
5036  basis%FACE_BASES(ni)%PTR=>new_sub_basis
5037  ENDIF
5038  ENDDO !ni
5039  ELSE
5040  ALLOCATE(basis%FACE_BASES(1),stat=err)
5041  IF(err/=0) CALL flagerror("Could not allocate basis face bases.",err,error,*999)
5042  basis%FACE_BASES(1)%PTR=>basis
5043  ENDIF
5044  ELSE
5045  ALLOCATE(basis%LINE_BASES(1),stat=err)
5046  IF(err/=0) CALL flagerror("Could not allocate basis line bases.",err,error,*999)
5047  basis%LINE_BASES(1)%PTR=>basis
5048  NULLIFY(basis%FACE_BASES)
5049  ENDIF
5050  ELSE
5051  CALL flagerror("Basis is not associated.",err,error,*999)
5052  ENDIF
5053 
5054  exits("BASIS_SIMPLEX_FAMILY_CREATE")
5055  RETURN
5056 999 IF(ASSOCIATED(new_sub_basis)) CALL basis_family_destroy(new_sub_basis%USER_NUMBER,new_sub_basis%FAMILY_NUMBER, &
5057  & dummy_err,dummy_error,*998)
5058 998 errorsexits("BASIS_SIMPLEX_FAMILY_CREATE",err,error)
5059  RETURN 1
5060  END SUBROUTINE basis_simplex_family_create
5061 
5062  !
5063  !================================================================================================================================
5064  !
5065 
5112  FUNCTION basis_simplex_basis_evaluate(BASIS,NODE_NUMBER,PARTIAL_DERIV_INDEX,XL,ERR,ERROR)
5113 
5114  !Argument variables
5115  TYPE(basis_type), POINTER :: BASIS
5116  INTEGER(INTG), INTENT(IN) :: NODE_NUMBER
5117  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIV_INDEX
5118  REAL(DP), INTENT(IN) :: XL(:)
5119  INTEGER(INTG), INTENT(OUT) :: ERR
5120  TYPE(varying_string), INTENT(OUT) :: ERROR
5121  !Function variable
5122  REAL(DP) :: BASIS_SIMPLEX_BASIS_EVALUATE
5123  !Local variables
5124  TYPE(varying_string) :: LOCAL_ERROR
5125 
5126  enters("BASIS_SIMPLEX_BASIS_EVALUATE",err,error,*999)
5127 
5128  basis_simplex_basis_evaluate=0.0_dp
5129  IF(ASSOCIATED(basis)) THEN
5130  IF(basis%TYPE==basis_simplex_type) THEN
5131  SELECT CASE(basis%NUMBER_OF_XI)
5132  CASE(1)
5133  SELECT CASE(partial_deriv_index)
5134  CASE(no_part_deriv)
5135  basis_simplex_basis_evaluate= &
5136  & basis_simplex_basis_derivative_evaluate(basis,node_number,no_part_deriv,xl,err,error)
5137  CASE(part_deriv_s1)
5138  basis_simplex_basis_evaluate= &
5139  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2,xl,err,error)
5140  IF(err/=0) GOTO 999
5141  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5142  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1,xl,err,error)
5143  CASE(part_deriv_s1_s1)
5144  basis_simplex_basis_evaluate= &
5145  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s1,xl,err,error)
5146  IF(err/=0) GOTO 999
5147  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5148  & 2.0_dp*basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s2,xl,err,error)
5149  IF(err/=0) GOTO 999
5150  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5151  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s2,xl,err,error)
5152  CASE DEFAULT
5153  local_error="The specified partial derivative index of "//trim(number_to_vstring(partial_deriv_index,"*",err,error))// &
5154  & " is invalid for a Simplex line basis."
5155  CALL flagerror(local_error,err,error,*999)
5156  END SELECT
5157  CASE(2)
5158  SELECT CASE(partial_deriv_index)
5159  CASE(no_part_deriv)
5160  basis_simplex_basis_evaluate= &
5161  & basis_simplex_basis_derivative_evaluate(basis,node_number,no_part_deriv,xl,err,error)
5162  CASE(part_deriv_s1)
5163  basis_simplex_basis_evaluate= &
5164  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3,xl,err,error)
5165  IF(err/=0) GOTO 999
5166  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5167  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1,xl,err,error)
5168  CASE(part_deriv_s1_s1)
5169  basis_simplex_basis_evaluate= &
5170  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s1,xl,err,error)
5171  IF(err/=0) GOTO 999
5172  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5173  & 2.0_dp*basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s3,xl,err,error)
5174  IF(err/=0) GOTO 999
5175  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5176  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s3,xl,err,error)
5177  CASE(part_deriv_s2)
5178  basis_simplex_basis_evaluate= &
5179  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3,xl,err,error)
5180  IF(err/=0) GOTO 999
5181  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5182  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2,xl,err,error)
5183  CASE(part_deriv_s2_s2)
5184  basis_simplex_basis_evaluate= &
5185  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s2,xl,err,error)
5186  IF(err/=0) GOTO 999
5187  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5188  & 2.0_dp*basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s3,xl,err,error)
5189  IF(err/=0) GOTO 999
5190  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5191  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s3,xl,err,error)
5192  CASE(part_deriv_s1_s2)
5193  basis_simplex_basis_evaluate= &
5194  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s3,xl,err,error)
5195  IF(err/=0) GOTO 999
5196  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5197  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s3,xl,err,error)
5198  IF(err/=0) GOTO 999
5199  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5200  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s3,xl,err,error)
5201  IF(err/=0) GOTO 999
5202  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5203  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s2,xl,err,error)
5204  CASE DEFAULT
5205  local_error="The specified partial derivative index of "//trim(number_to_vstring(partial_deriv_index,"*",err,error))// &
5206  & " is invalid for a Simplex triangle basis."
5207  CALL flagerror(local_error,err,error,*999)
5208  END SELECT
5209  CASE(3)
5210  SELECT CASE(partial_deriv_index)
5211  CASE(no_part_deriv)
5212  basis_simplex_basis_evaluate= &
5213  & basis_simplex_basis_derivative_evaluate(basis,node_number,no_part_deriv,xl,err,error)
5214  CASE(part_deriv_s1)
5215  basis_simplex_basis_evaluate= &
5216  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4,xl,err,error)
5217  IF(err/=0) GOTO 999
5218  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5219  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1,xl,err,error)
5220  CASE(part_deriv_s1_s1)
5221  basis_simplex_basis_evaluate= &
5222  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s1,xl,err,error)
5223  IF(err/=0) GOTO 999
5224  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5225  & 2.0_dp*basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s4,xl,err,error)
5226  IF(err/=0) GOTO 999
5227  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5228  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4_s4,xl,err,error)
5229  CASE(part_deriv_s2)
5230  basis_simplex_basis_evaluate= &
5231  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4,xl,err,error)
5232  IF(err/=0) GOTO 999
5233  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5234  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2,xl,err,error)
5235  CASE(part_deriv_s2_s2)
5236  basis_simplex_basis_evaluate= &
5237  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s2,xl,err,error)
5238  IF(err/=0) GOTO 999
5239  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5240  & 2.0_dp*basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s4,xl,err,error)
5241  IF(err/=0) GOTO 999
5242  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5243  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4_s4,xl,err,error)
5244  CASE(part_deriv_s1_s2)
5245  basis_simplex_basis_evaluate= &
5246  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4_s4,xl,err,error)
5247  IF(err/=0) GOTO 999
5248  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5249  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s4,xl,err,error)
5250  IF(err/=0) GOTO 999
5251  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5252  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s4,xl,err,error)
5253  IF(err/=0) GOTO 999
5254  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5255  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s2,xl,err,error)
5256  CASE(part_deriv_s3)
5257  basis_simplex_basis_evaluate= &
5258  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4,xl,err,error)
5259  IF(err/=0) GOTO 999
5260  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5261  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3,xl,err,error)
5262  CASE(part_deriv_s3_s3)
5263  basis_simplex_basis_evaluate= &
5264  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s3,xl,err,error)
5265  IF(err/=0) GOTO 999
5266  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5267  & 2.0_dp*basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s4,xl,err,error)
5268  IF(err/=0) GOTO 999
5269  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5270  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4_s4,xl,err,error)
5271  CASE(part_deriv_s1_s3)
5272  basis_simplex_basis_evaluate= &
5273  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4_s4,xl,err,error)
5274  IF(err/=0) GOTO 999
5275  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5276  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s4,xl,err,error)
5277  IF(err/=0) GOTO 999
5278  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5279  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s4,xl,err,error)
5280  IF(err/=0) GOTO 999
5281  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5282  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s3,xl,err,error)
5283  CASE(part_deriv_s2_s3)
5284  basis_simplex_basis_evaluate= &
5285  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4_s4,xl,err,error)
5286  IF(err/=0) GOTO 999
5287  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5288  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s4,xl,err,error)
5289  IF(err/=0) GOTO 999
5290  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5291  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s4,xl,err,error)
5292  IF(err/=0) GOTO 999
5293  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5294  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s3,xl,err,error)
5295  CASE(part_deriv_s1_s2_s3)
5296  basis_simplex_basis_evaluate= &
5297  basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s4_s4_s4,xl,err,error)
5298  IF(err/=0) GOTO 999
5299  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5300  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s4_s4,xl,err,error)
5301  IF(err/=0) GOTO 999
5302  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5303  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s4_s4,xl,err,error)
5304  IF(err/=0) GOTO 999
5305  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5306  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s3_s4_s4,xl,err,error)
5307  IF(err/=0) GOTO 999
5308  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5309  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s2_s4,xl,err,error)
5310  IF(err/=0) GOTO 999
5311  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5312  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s3_s4,xl,err,error)
5313  IF(err/=0) GOTO 999
5314  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate+ &
5315  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s2_s3_s4,xl,err,error)
5316  IF(err/=0) GOTO 999
5317  basis_simplex_basis_evaluate=basis_simplex_basis_evaluate- &
5318  & basis_simplex_basis_derivative_evaluate(basis,node_number,part_deriv_s1_s2_s3,xl,err,error)
5319  CASE DEFAULT
5320  local_error="The specified partial derivative index of "//trim(number_to_vstring(partial_deriv_index,"*",err,error))// &
5321  & " is invalid for a Simplex tetrahedra basis."
5322  CALL flagerror(local_error,err,error,*999)
5323  END SELECT
5324  CASE DEFAULT
5325  local_error="Invalid number of Xi coordinates. The number of xi coordinates for this basis is "// &
5326  & trim(number_to_vstring(basis%NUMBER_OF_XI,"*",err,error))//" which should be between 1 and 3."
5327  CALL flagerror(local_error,err,error,*999)
5328  END SELECT
5329  IF(err/=0) GOTO 999
5330  ELSE
5331  CALL flagerror("Basis is not a simplex basis.",err,error,*999)
5332  ENDIF
5333  ELSE
5334  CALL flagerror("Basis is not associated.",err,error,*999)
5335  ENDIF
5336 
5337  exits("BASIS_SIMPLEX_BASIS_EVALUATE.")
5338  RETURN
5339 999 errorsexits("BASIS_SIMPLEX_BASIS_EVALUATE.",err,error)
5340  RETURN
5341  END FUNCTION basis_simplex_basis_evaluate
5342 
5343  !
5344  !================================================================================================================================
5345  !
5346 
5348  FUNCTION basis_simplex_basis_derivative_evaluate(BASIS,NODE_NUMBER,PARTIAL_DERIV_INDEX,XL,ERR,ERROR)
5349 
5350  !Argument variables
5351  TYPE(basis_type), POINTER :: BASIS
5352  INTEGER(INTG), INTENT(IN) :: NODE_NUMBER
5353  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIV_INDEX
5354  REAL(DP), INTENT(IN) :: XL(:)
5355  INTEGER(INTG), INTENT(OUT) :: ERR
5356  TYPE(varying_string), INTENT(OUT) :: ERROR
5357  !Function variable
5358  REAL(DP) :: BASIS_SIMPLEX_BASIS_DERIVATIVE_EVALUATE
5359  !Local variables
5360  INTEGER(INTG) :: nic
5361  TYPE(varying_string) :: LOCAL_ERROR
5362 
5363  enters("BASIS_SIMPLEX_BASIS_DERIVATIVE_EVALUATE",err,error,*999)
5364 
5365  basis_simplex_basis_derivative_evaluate=1.0_dp
5366  IF(ASSOCIATED(basis)) THEN
5367  DO nic=1,basis%NUMBER_OF_XI_COORDINATES
5368  SELECT CASE(basis%INTERPOLATION_ORDER(nic))
5369  CASE(basis_linear_interpolation_order)
5370  basis_simplex_basis_derivative_evaluate=basis_simplex_basis_derivative_evaluate* &
5371  & simplex_linear_evaluate(basis%NODE_POSITION_INDEX(node_number,nic), &
5372  & partial_derivative_index(partial_deriv_index,nic),xl(nic),err,error)
5373  CASE(basis_quadratic_interpolation_order)
5374  basis_simplex_basis_derivative_evaluate=basis_simplex_basis_derivative_evaluate* &
5375  & simplex_quadratic_evaluate(basis%NODE_POSITION_INDEX(node_number,nic), &
5376  & partial_derivative_index(partial_deriv_index,nic),xl(nic),err,error)
5377  CASE(basis_cubic_interpolation_order)
5378  basis_simplex_basis_derivative_evaluate=basis_simplex_basis_derivative_evaluate* &
5379  & simplex_cubic_evaluate(basis%NODE_POSITION_INDEX(node_number,nic), &
5380  & partial_derivative_index(partial_deriv_index,nic),xl(nic),err,error)
5381  CASE DEFAULT
5382  local_error="Interpolation order value "//trim(number_to_vstring(basis%INTERPOLATION_ORDER(nic),"*",err,error))// &
5383  & " for xi coordinate direction "//trim(number_to_vstring(nic,"*",err,error))//" is invalid."
5384  CALL flagerror(local_error,err,error,*999)
5385  END SELECT
5386  IF(err/=0) GOTO 999
5387  ENDDO !nic
5388  ELSE
5389  CALL flagerror("Basis is not associated.",err,error,*999)
5390  ENDIF
5391 
5392  exits("BASIS_SIMPLEX_BASIS_DERIVATIVE_EVALUATE")
5393  RETURN
5394 999 errorsexits("BASIS_SIMPLEX_BASIS_DERIVATIVE_EVALUATE",err,error)
5395  RETURN
5396  END FUNCTION basis_simplex_basis_derivative_evaluate
5397 
5398  !
5399  !================================================================================================================================
5400  !
5401 
5403  SUBROUTINE basis_sub_basis_create(PARENT_BASIS,NUMBER_OF_XI,XI_DIRECTIONS,SUB_BASIS,ERR,ERROR,*)
5404 
5405  !Argument variables
5406  TYPE(basis_type), POINTER :: PARENT_BASIS
5407  INTEGER(INTG), INTENT(IN) :: NUMBER_OF_XI
5408  INTEGER(INTG), INTENT(IN) :: XI_DIRECTIONS(:)
5409  TYPE(basis_type), POINTER :: SUB_BASIS
5410  INTEGER(INTG), INTENT(OUT) :: ERR
5411  TYPE(varying_string), INTENT(OUT) :: ERROR
5412  !Local Variables
5413  INTEGER(INTG) :: basis_idx,ni,NUMBER_COLLAPSED,NUMBER_END_COLLAPSED
5414  TYPE(basis_type), POINTER :: NEW_SUB_BASIS
5415  TYPE(basis_ptr_type), POINTER :: NEW_SUB_BASES(:)
5416  TYPE(varying_string) :: LOCAL_ERROR
5417 
5418  NULLIFY(new_sub_basis)
5419  NULLIFY(new_sub_bases)
5420 
5421  enters("BASIS_SUB_BASIS_CREATE",err,error,*999)
5422 
5423  IF(ASSOCIATED(parent_basis)) THEN
5424  IF(ASSOCIATED(sub_basis)) THEN
5425  CALL flagerror("The sub-basis is already associated",err,error,*999)
5426  ELSE
5427  IF(number_of_xi>0.AND.number_of_xi<4) THEN
5428  IF(any(xi_directions<1).OR.any(xi_directions>3)) CALL flagerror("Invalid xi directions specified",err,error,*999)
5429  IF(SIZE(xi_directions,1)/=number_of_xi) &
5430  & CALL flagerror("The size of the xi directions array must be the same as the number of xi directions",err,error,*999)
5431  ALLOCATE(new_sub_basis,stat=err)
5432  IF(err/=0) CALL flagerror("Could not allocate sub-basis",err,error,*999)
5433  new_sub_basis%USER_NUMBER=parent_basis%USER_NUMBER
5434  new_sub_basis%GLOBAL_NUMBER=parent_basis%GLOBAL_NUMBER
5435  new_sub_basis%FAMILY_NUMBER=parent_basis%NUMBER_OF_SUB_BASES+1
5436  new_sub_basis%NUMBER_OF_SUB_BASES=0
5437  NULLIFY(new_sub_basis%SUB_BASES)
5438  new_sub_basis%PARENT_BASIS=>parent_basis
5439  new_sub_basis%NUMBER_OF_XI=number_of_xi
5440  new_sub_basis%TYPE=parent_basis%TYPE
5441  ALLOCATE(new_sub_basis%INTERPOLATION_XI(number_of_xi),stat=err)
5442  IF(err/=0) CALL flagerror("Could not allocate sub-basis interpolation xi",err,error,*999)
5443  ALLOCATE(new_sub_basis%COLLAPSED_XI(number_of_xi),stat=err)
5444  IF(err/=0) CALL flagerror("Could not allocate sub-basis collapsed xi",err,error,*999)
5445  number_collapsed=0
5446  number_end_collapsed=0
5447  DO ni=1,number_of_xi
5448  new_sub_basis%INTERPOLATION_XI(ni)=parent_basis%INTERPOLATION_XI(xi_directions(ni))
5449  new_sub_basis%COLLAPSED_XI(ni)=parent_basis%COLLAPSED_XI(xi_directions(ni))
5450  IF(new_sub_basis%COLLAPSED_XI(ni)==basis_xi_collapsed) THEN
5451  number_collapsed=number_collapsed+1
5452  ELSE IF(new_sub_basis%COLLAPSED_XI(ni)==basis_collapsed_at_xi0.OR.new_sub_basis%COLLAPSED_XI(ni)== &
5453  & basis_collapsed_at_xi1) THEN
5454  number_end_collapsed=number_end_collapsed+1
5455  ENDIF
5456  ENDDO !ni
5457  IF(number_collapsed==0.OR.number_end_collapsed==0) new_sub_basis%COLLAPSED_XI(1:number_of_xi)=basis_not_collapsed
5458  NULLIFY(new_sub_basis%QUADRATURE%BASIS)
5459  CALL basis_quadrature_initialise(new_sub_basis,err,error,*999)
5460  new_sub_basis%QUADRATURE%TYPE=parent_basis%QUADRATURE%TYPE
5461  DO ni=1,number_of_xi
5462  new_sub_basis%QUADRATURE%NUMBER_OF_GAUSS_XI(ni)=parent_basis%QUADRATURE%NUMBER_OF_GAUSS_XI(xi_directions(ni))
5463  ENDDO !ni
5464  new_sub_basis%BASIS_FINISHED=.true.
5465  IF(number_of_xi>1) THEN
5466  ALLOCATE(new_sub_basis%LINE_BASES(number_of_xi),stat=err)
5467  IF(err/=0) CALL flagerror("Could not allocate sub-basis line bases",err,error,*999)
5468  IF(number_of_xi>2) THEN
5469  ALLOCATE(new_sub_basis%FACE_BASES(number_of_xi),stat=err)
5470  IF(err/=0) CALL flagerror("Could not allocate sub-basis face bases",err,error,*999)
5471  ELSE
5472  ALLOCATE(new_sub_basis%FACE_BASES(1),stat=err)
5473  IF(err/=0) CALL flagerror("Could not allocate sub-basis face bases",err,error,*999)
5474  new_sub_basis%FACE_BASES(1)%PTR=>new_sub_basis
5475  ENDIF
5476  ELSE
5477  ALLOCATE(new_sub_basis%LINE_BASES(1),stat=err)
5478  IF(err/=0) CALL flagerror("Could not allocate basis line bases",err,error,*999)
5479  new_sub_basis%LINE_BASES(1)%PTR=>new_sub_basis
5480  NULLIFY(new_sub_basis%FACE_BASES)
5481  ENDIF
5482  !Add the new sub-basis to the list of sub-bases in the parent basis
5483  ALLOCATE(new_sub_bases(parent_basis%NUMBER_OF_SUB_BASES+1),stat=err)
5484  IF(err/=0) CALL flagerror("Could not allocate new sub-bases",err,error,*999)
5485  DO basis_idx=1,parent_basis%NUMBER_OF_SUB_BASES
5486  new_sub_bases(basis_idx)%PTR=>parent_basis%SUB_BASES(basis_idx)%PTR
5487  ENDDO !basis_idx
5488  new_sub_bases(parent_basis%NUMBER_OF_SUB_BASES+1)%PTR=>new_sub_basis
5489  parent_basis%NUMBER_OF_SUB_BASES=parent_basis%NUMBER_OF_SUB_BASES+1
5490  IF(ASSOCIATED(parent_basis%SUB_BASES)) DEALLOCATE(parent_basis%SUB_BASES)
5491  parent_basis%SUB_BASES=>new_sub_bases
5492  sub_basis=>new_sub_basis
5493  ELSE
5494  local_error="Invalid number of xi directions specified ("// &
5495  & trim(number_to_vstring(number_of_xi,"*",err,error))//"). You must specify between 1 and 3 xi directions"
5496  CALL flagerror(local_error,err,error,*999)
5497  ENDIF
5498  ENDIF
5499  ELSE
5500  CALL flagerror("Parent basis is not associated",err,error,*999)
5501  ENDIF
5502 
5503  exits("BASIS_SUB_BASIS_CREATE")
5504  RETURN
5505 999 errorsexits("BASIS_SUB_BASIS_CREATE",err,error)
5506  RETURN 1
5507  END SUBROUTINE basis_sub_basis_create
5508 
5509  !
5510  !================================================================================================================================
5511  !
5512 
5514  SUBROUTINE basis_type_get(BASIS,TYPE,ERR,ERROR,*)
5515 
5516  !Argument variables
5517  TYPE(basis_type), POINTER :: BASIS
5518  INTEGER(INTG), INTENT(OUT) :: TYPE
5519  INTEGER(INTG), INTENT(OUT) :: ERR
5520  TYPE(varying_string), INTENT(OUT) :: ERROR
5521  !Local Variables
5522 
5523  enters("BASIS_TYPE_GET",err,error,*999)
5524 
5525  IF(ASSOCIATED(basis)) THEN
5526  IF(basis%BASIS_FINISHED) THEN
5527  TYPE=basis%TYPE
5528  ELSE
5529  CALL flagerror("Basis has not been finished yet",err,error,*999)
5530  ENDIF
5531  ELSE
5532  CALL flagerror("Basis is not associated",err,error,*999)
5533  ENDIF
5534 
5535  exits("BASIS_TYPE_GET")
5536  RETURN
5537 999 errorsexits("BASIS_TYPE_GET",err,error)
5538  RETURN
5539  END SUBROUTINE basis_type_get
5540 
5541  !
5542  !================================================================================================================================
5543  !
5544 
5546  SUBROUTINE basis_type_set_number(USER_NUMBER,TYPE,ERR,ERROR,*)
5547 
5548  !Argument variables
5549  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
5550  INTEGER(INTG), INTENT(IN) :: TYPE
5551  INTEGER(INTG), INTENT(OUT) :: ERR
5552  TYPE(varying_string), INTENT(OUT) :: ERROR
5553  !Local Variables
5554  TYPE(basis_type), POINTER :: BASIS
5555 
5556  enters("BASIS_TYPE_SET_NUMBER",err,error,*999)
5557 
5558  CALL basis_user_number_find(user_number,basis,err,error,*999)
5559  CALL basis_type_set_ptr(basis,TYPE,ERR,ERROR,*999)
5560 
5561  exits("BASIS_TYPE_SET_NUMBER")
5562  RETURN
5563 999 errorsexits("BASIS_TYPE_SET_NUMBER",err,error)
5564  RETURN 1
5565  END SUBROUTINE basis_type_set_number
5566 
5567  !
5568  !================================================================================================================================
5569  !
5570 
5572  SUBROUTINE