Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_interpolation_hex_hermite.F90
Go to the documentation of this file.
1 module sll_m_interpolation_hex_hermite
2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 #include "sll_memory.h"
4 #include "sll_working_precision.h"
5 
6  use sll_m_hermite_interpolation_2d, only: &
7  sll_s_compute_w_hermite
8 
9  use sll_m_hexagonal_meshes, only: &
10  sll_s_get_cell_vertices_index, &
11  sll_s_get_edge_index, &
12  sll_s_get_triangle_index, &
13  sll_t_hex_mesh_2d
14 
15  implicit none
16 
17  public :: &
18  sll_s_der_finite_difference, &
19  sll_s_hermite_interpolation, &
20  sll_s_print_method
21 
22  private
23 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
24 
25 contains
26 
27  subroutine sll_s_der_finite_difference(f_tn, p, step, mesh, deriv)
28  !-> computation of the partial derivatives in the directions H1, H2 and H3
29  ! with dirichlet boundary condition
30  type(sll_t_hex_mesh_2d), pointer :: mesh
31  sll_real64, dimension(:, :), intent(out):: deriv
32  sll_real64, dimension(:), intent(in) :: f_tn
33  sll_real64, intent(in) :: step
34  sll_int32, intent(in) :: p
35  sll_int32 :: num_cells, i, j, r, s, n_points
36  sll_int32 :: n1, n2, n3, n4, n5, n6, nc1
37  sll_real64 :: dh1, dh2, dh3, dh4, dh5, dh6
38  sll_int32 :: h1, h2, h1t, h2t, p2
39  sll_int32 :: index_boundary
40  sll_int32 :: index_boundary1
41  sll_int32 :: index_boundary2
42  sll_int32 :: index_boundary3
43  sll_int32 :: index_boundary4
44  sll_int32 :: index_boundary5
45  sll_int32 :: index_boundary6
46  sll_int32 :: maxrs
47  logical :: opsign1
48  logical :: opsign2
49  logical :: opsign3
50  logical :: boundary
51  logical :: boundary1
52  logical :: boundary2
53  logical :: boundary3
54  logical :: boundary4
55  logical :: boundary5
56  logical :: boundary6
57  logical :: secure
58  sll_real64, dimension(:), allocatable :: w
59  sll_real64, dimension(:, :), allocatable :: f
60 
61  n_points = size(f_tn)
62  num_cells = mesh%num_cells
63 
64  !*********************************************************************
65  ! computation of the coefficients in fonction of the degree p
66  ! of the approximation required
67  !*********************************************************************
68 
69  p2 = p/2
70 
71  if ((p2)*2 == p) then !if p is even
72  r = -p2
73  s = p2
74  else !if p is odd
75  r = -p2 - 1 ! de-centered to the left
76  s = p2 !+ 1 ! de-centered to the right
77  end if
78 
79  allocate (w(r:s), f(1:6, r:s))
80 
81  call sll_s_compute_w_hermite(w, r, s)
82 
83  !***********************************************************************
84  ! Computation of the derivatives in both hexa directions h1 and h2,
85  ! then for both x and y directions for every points of the mesh
86  !***********************************************************************
87 
88  maxrs = max(abs(r), s)
89  nc1 = num_cells + 1
90 
91  do i = 1, n_points
92 
93  secure = .false.
94  h1 = mesh%hex_coord(1, i)
95  h2 = mesh%hex_coord(2, i)
96 
97  ! checking if the point "i" is in the "secure zone"
98  ! where boundary conditions won't play a role
99 
100  if (h1 <= num_cells - maxrs .and. h2 <= num_cells - maxrs &
101  .and. h1 >= -num_cells + maxrs .and. h2 >= -num_cells + maxrs) then
102  if (h1 <= 0 .and. h2 >= 0 .and. h2 - h1 <= num_cells - maxrs &
103  .or. h2 <= 0 .and. h1 >= 0 .and. h1 - h2 <= num_cells - maxrs &
104  ) then
105  secure = .true.
106  end if
107  end if
108 
109  if (secure) then
110 
111  do j = r, s !find the required points in the h1, h2 and h3 directions
112 
113  h1t = h1 + j
114  h2t = h2 + j
115 
116  n1 = mesh%hex_to_global(h1t, h2)
117  f(1, j) = f_tn(n1)
118 
119  n2 = mesh%hex_to_global(h1, h2t)
120  f(2, j) = f_tn(n2)
121 
122  n3 = mesh%hex_to_global(h1t, h2t)
123  f(3, j) = f_tn(n3)
124 
125  !find the required points in the -h1, -h2 and -h3 directions
126 
127  h1t = h1 - j
128  h2t = h2 - j
129 
130  n4 = mesh%hex_to_global(h1t, h2)
131  f(4, j) = f_tn(n4)
132 
133  n5 = mesh%hex_to_global(h1, h2t)
134  f(5, j) = f_tn(n5)
135 
136  n6 = mesh%hex_to_global(h1t, h2t)
137  f(6, j) = f_tn(n6)
138 
139  end do
140 
141  else
142 
143  boundary = .false.
144 
145  ! checking if the point "i" is on the boundary:
146 
147  if (h1 == num_cells .and. h2 >= 0 .and. h2 <= num_cells & !edge 1
148  .or. h1 == -num_cells .and. h2 <= 0 .and. h2 >= -num_cells & !edge4
149  .or. h2 == num_cells .and. h1 >= 0 .and. h1 <= num_cells & !edge2
150  .or. h2 == -num_cells .and. h1 <= 0 .and. h1 >= -num_cells & !edge5
151  .or. h1 <= 0 .and. h2 >= 0 .and. h2 - h1 == num_cells & !edge3
152  .or. h2 <= 0 .and. h1 >= 0 .and. h1 - h2 == num_cells & !edge6
153  ) then
154  boundary = .true.
155  index_boundary = i
156  end if
157 
158  ! we could have kept in memory which edge it is on, but due to simplicity concern
159  ! we did not implement a test for each edge even though it would save tests
160  ! it could be done for optimisation purpose in the future
161 
162  if (boundary) then
163 
164  do j = r, s !find the required points in the h1, h2 and h3 directions
165 
166  opsign1 = .false.
167  opsign2 = .false.
168  opsign3 = .false.
169 
170  h1t = h1 + j
171  h2t = h2 + j
172 
173  if (h2*h1t < 0) opsign1 = .true.
174  if (h2t*h1 < 0) opsign2 = .true.
175  if (h2t*h1t < 0) opsign3 = .true.
176 
177  !test if outside mesh
178 
179  if (abs(h1t) > num_cells .or. opsign1 .and. (abs(h1t) + abs(h2) > num_cells)) then
180  f(1, j) = f_tn(index_boundary) ! dirichlet boundary condition
181  else
182  n1 = mesh%hex_to_global(h1t, h2)
183  f(1, j) = f_tn(n1)
184  end if
185 
186  if (abs(h2t) > num_cells .or. opsign2 .and. (abs(h1) + abs(h2t) > num_cells)) then
187  f(2, j) = f_tn(index_boundary)
188  else
189  n2 = mesh%hex_to_global(h1, h2t)
190  f(2, j) = f_tn(n2)
191  end if
192 
193  if (abs(h1t) > num_cells .or. abs(h2t) > num_cells .or. opsign3 .and. (abs(h1t) + abs(h2t) > num_cells)) then
194  f(3, j) = f_tn(index_boundary)
195  else
196  n3 = mesh%hex_to_global(h1t, h2t)
197  f(3, j) = f_tn(n3)
198  end if
199 
200  !find the required points in the -h1, -h2 and -h3 directions
201 
202  opsign1 = .false.
203  opsign2 = .false.
204  opsign3 = .false.
205 
206  h1t = h1 - j
207  h2t = h2 - j
208 
209  if (h2*h1t < 0) opsign1 = .true.
210  if (h2t*h1 < 0) opsign2 = .true.
211  if (h2t*h1t < 0) opsign3 = .true.
212 
213  !test if outside mesh
214 
215  if (abs(h1t) > num_cells .or. &
216  opsign1 .and. (abs(h1t) + abs(h2) > num_cells)) then
217  f(4, j) = f_tn(index_boundary) ! dirichlet boundary condition
218  else
219  n4 = mesh%hex_to_global(h1t, h2)
220  f(4, j) = f_tn(n4)
221  end if
222 
223  if (abs(h2t) > num_cells .or. &
224  opsign2 .and. (abs(h1) + abs(h2t) > num_cells)) then
225  f(5, j) = f_tn(index_boundary)
226  else
227  n5 = mesh%hex_to_global(h1, h2t)
228  f(5, j) = f_tn(n5)
229  end if
230 
231  if (abs(h1t) > num_cells .or. abs(h2t) > num_cells .or. &
232  opsign3 .and. (abs(h1t) + abs(h2t) > num_cells)) then
233  f(6, j) = f_tn(index_boundary)
234  else
235  n6 = mesh%hex_to_global(h1t, h2t)
236  f(6, j) = f_tn(n6)
237  end if
238 
239  end do
240 
241  else !case where there will be boundary points but "i" is not one of them
242 
243  f(1:6, 0) = f_tn(i)
244 
245  boundary1 = .false.
246  boundary2 = .false.
247  boundary3 = .false.
248  boundary4 = .false.
249  boundary5 = .false.
250  boundary6 = .false.
251 
252  do j = 1, s
253 
254  h1t = h1 + j
255  h2t = h2 + j
256 
257  ! we need to stock the index of the boundary (ies) point(s) in any direction
258  ! then test
259 
260  if (boundary1 .eqv. .false.) then
261  if (h1t == num_cells .and. h2 >= 0 .and. h2 <= num_cells .or. &
262  h2 <= 0 .and. h1t >= 0 .and. h1t - h2 == num_cells) then
263  boundary1 = .true.
264  index_boundary1 = mesh%hex_to_global(h1t, h2)
265  end if
266  end if
267 
268  if (boundary2 .eqv. .false.) then
269  if (h2t == num_cells .and. h1 >= 0 .and. h1 <= num_cells .or. &
270  h1 <= 0 .and. h2t >= 0 .and. h2t - h1 == num_cells) then
271  boundary2 = .true.
272  index_boundary2 = mesh%hex_to_global(h1, h2t)
273  end if
274  end if
275 
276  if (boundary3 .eqv. .false.) then
277  if (h1t == num_cells .and. h2t >= 0 .and. h2t <= num_cells .or. &
278  h2t == num_cells .and. h1t >= 0 .and. h1t <= num_cells) then
279  boundary3 = .true.
280  index_boundary3 = mesh%hex_to_global(h1t, h2t)
281  end if
282  end if
283 
284  if (boundary1) then
285  f(1, j) = f_tn(index_boundary1) ! dirichlet boundary condition
286  else
287  n1 = mesh%hex_to_global(h1t, h2)
288  f(1, j) = f_tn(n1)
289  end if
290 
291  if (boundary2) then
292  f(2, j) = f_tn(index_boundary2)
293  else
294  n2 = mesh%hex_to_global(h1, h2t)
295  f(2, j) = f_tn(n2)
296  end if
297 
298  if (boundary3) then
299  f(3, j) = f_tn(index_boundary3)
300  else
301  n3 = mesh%hex_to_global(h1t, h2t)
302  f(3, j) = f_tn(n3)
303  end if
304 
305  h1t = h1 - j
306  h2t = h2 - j
307 
308  if (boundary4 .eqv. .false.) then
309  if (h1t == -num_cells .and. h2 <= 0 .and. h2 >= -num_cells .or. &
310  h1t <= 0 .and. h2 >= 0 .and. h2 - h1t == num_cells) then
311  boundary4 = .true.
312  index_boundary4 = mesh%hex_to_global(h1t, h2)
313  end if
314  end if
315 
316  if (boundary5 .eqv. .false.) then
317  if (h2t == -num_cells .and. h1 <= 0 .and. h1 >= -num_cells .or. &
318  h2t <= 0 .and. h1 >= 0 .and. h1 - h2t == num_cells) then
319  boundary5 = .true.
320  index_boundary5 = mesh%hex_to_global(h1, h2t)
321  end if
322  end if
323 
324  if (boundary6 .eqv. .false.) then
325  if (h1t == -num_cells .and. h2t <= 0 .and. h2t >= -num_cells .or. &
326  h2t == -num_cells .and. h1t <= 0 .and. h1t >= -num_cells) then
327  boundary6 = .true.
328  index_boundary6 = mesh%hex_to_global(h1t, h2t)
329  end if
330  end if
331 
332  if (boundary4) then
333  f(4, j) = f_tn(index_boundary4) ! dirichlet boundary condition
334  else
335  n4 = mesh%hex_to_global(h1t, h2)
336  f(4, j) = f_tn(n4)
337  end if
338 
339  if (boundary5) then
340  f(5, j) = f_tn(index_boundary5)
341  else
342  n5 = mesh%hex_to_global(h1, h2t)
343  f(5, j) = f_tn(n5)
344  end if
345 
346  if (boundary6) then
347  f(6, j) = f_tn(index_boundary6)
348  else
349  n6 = mesh%hex_to_global(h1t, h2t)
350  f(6, j) = f_tn(n6)
351  end if
352 
353  end do
354 
355  boundary1 = .false.
356  boundary2 = .false.
357  boundary3 = .false.
358  boundary4 = .false.
359  boundary5 = .false.
360  boundary6 = .false.
361 
362  do j = -1, r, -1
363 
364  h1t = h1 + j
365  h2t = h2 + j
366 
367  if (boundary1 .eqv. .false.) then
368  if (h1t == -num_cells .and. h2 <= 0 .and. h2 >= -num_cells .or. &
369  h1t <= 0 .and. h2 >= 0 .and. h2 - h1t == num_cells) then
370  boundary1 = .true.
371  index_boundary1 = mesh%hex_to_global(h1t, h2)
372  end if
373  end if
374 
375  if (boundary2 .eqv. .false.) then
376  if (h2t == -num_cells .and. h1 <= 0 .and. h1 >= -num_cells .or. &
377  h2t <= 0 .and. h1 >= 0 .and. h1 - h2t == num_cells) then
378  boundary2 = .true.
379  index_boundary2 = mesh%hex_to_global(h1, h2t)
380  end if
381  end if
382 
383  if (boundary3 .eqv. .false.) then
384  if (h1t == -num_cells .and. h2t <= 0 .and. h2t >= -num_cells .or. &
385  h2t == -num_cells .and. h1t <= 0 .and. h1t >= -num_cells) then
386  boundary3 = .true.
387  index_boundary3 = mesh%hex_to_global(h1t, h2t)
388  end if
389  end if
390 
391  if (boundary1) then
392  f(1, j) = f_tn(index_boundary1) ! dirichlet boundary condition
393  else
394  n1 = mesh%hex_to_global(h1t, h2)
395  f(1, j) = f_tn(n1)
396  end if
397 
398  if (boundary2) then
399  f(2, j) = f_tn(index_boundary2)
400  else
401  n2 = mesh%hex_to_global(h1, h2t)
402  f(2, j) = f_tn(n2)
403  end if
404 
405  if (boundary3) then
406  f(3, j) = f_tn(index_boundary3)
407  else
408  n3 = mesh%hex_to_global(h1t, h2t)
409  f(3, j) = f_tn(n3)
410  end if
411 
412  h1t = h1 - j
413  h2t = h2 - j
414 
415  if (boundary4 .eqv. .false.) then
416  if (h1t == num_cells .and. h2 >= 0 .and. h2 <= num_cells .or. &
417  h2 <= 0 .and. h1t >= 0 .and. h1t - h2 == num_cells) then
418  boundary4 = .true.
419  index_boundary4 = mesh%hex_to_global(h1t, h2)
420  end if
421  end if
422 
423  if (boundary5 .eqv. .false.) then
424  if (h2t == num_cells .and. h1 >= 0 .and. h1 <= num_cells .or. &
425  h1 <= 0 .and. h2t >= 0 .and. h2t - h1 == num_cells) then
426  boundary5 = .true.
427  index_boundary5 = mesh%hex_to_global(h1, h2t)
428  end if
429  end if
430 
431  if (boundary6 .eqv. .false.) then
432  if (h1t == num_cells .and. h2t >= 0 .and. h2t <= num_cells .or. &
433  h2t == num_cells .and. h1t >= 0 .and. h1t <= num_cells) then
434  boundary6 = .true.
435  index_boundary6 = mesh%hex_to_global(h1t, h2t)
436  end if
437  end if
438 
439  if (boundary4) then
440  f(4, j) = f_tn(index_boundary4) ! dirichlet boundary condition
441  else
442  n4 = mesh%hex_to_global(h1t, h2)
443  f(4, j) = f_tn(n4)
444  end if
445 
446  if (boundary5) then
447  f(5, j) = f_tn(index_boundary5)
448  else
449  n5 = mesh%hex_to_global(h1, h2t)
450  f(5, j) = f_tn(n5)
451  end if
452 
453  if (boundary6) then
454  f(6, j) = f_tn(index_boundary6)
455  else
456  n6 = mesh%hex_to_global(h1t, h2t)
457  f(6, j) = f_tn(n6)
458  end if
459 
460  end do
461 
462  end if
463 
464  end if
465 
466  dh1 = 0._f64
467  dh2 = 0._f64
468  dh3 = 0._f64
469  dh4 = 0._f64
470  dh5 = 0._f64
471  dh6 = 0._f64
472 
473  do j = r, s
474  dh1 = dh1 + w(j)*f(1, j)
475  dh2 = dh2 + w(j)*f(2, j)
476  dh3 = dh3 + w(j)*f(3, j)
477  dh4 = dh4 + w(j)*f(4, j)
478  dh5 = dh5 + w(j)*f(5, j)
479  dh6 = dh6 + w(j)*f(6, j)
480  end do
481 
482  deriv(1, i) = dh1/step
483  deriv(2, i) = dh2/step
484  deriv(3, i) = dh3/step
485  deriv(4, i) = dh4/step
486  deriv(5, i) = dh5/step
487  deriv(6, i) = dh6/step
488 
489  end do
490 
491  deallocate (f, w)
492 
493  end subroutine sll_s_der_finite_difference
494 
495  subroutine sll_s_hermite_interpolation(num, x, y, f_tn, center_value, edge_value, output_tn1, mesh, deriv, aire, num_method)
496 
497  type(sll_t_hex_mesh_2d), pointer :: mesh
498  sll_real64, dimension(:), intent(in) :: f_tn
499  sll_real64, dimension(:), intent(in) :: edge_value
500  sll_real64, dimension(:), intent(in) :: center_value
501  sll_real64, dimension(:, :), intent(in) :: deriv
502  sll_real64, dimension(:), intent(out) :: output_tn1
503  sll_real64, intent(in) :: x, y, aire
504  sll_int32, intent(in) :: num, num_method
505  sll_real64 :: x1, x2, x3, y1, y2, y3, f, step
506  sll_real64 :: l1, l2, l3
507  sll_real64, dimension(:), allocatable :: freedom, base
508  sll_real64 :: a2
509  sll_real64 :: x1x, x2x, x3x, y1y, y2y, y3y
510  sll_int32 :: i, i1, i2, i3, k11, k12, center_index
511  sll_int32 :: num_degree
512  sll_int32 :: edge_index1, edge_index2, edge_index3
513  logical :: test
514 
515  if (num_method == 9) then
516  allocate (freedom(9), base(9))
517  num_degree = 9
518  elseif (num_method == 10) then
519  allocate (freedom(10), base(10))
520  num_degree = 10
521  elseif (num_method == 11) then
522  allocate (freedom(9), base(9))
523  num_degree = 9
524  elseif (num_method == 12) then
525  allocate (freedom(12), base(12))
526  num_degree = 12
527  elseif (num_method == 15) then
528  allocate (freedom(15), base(15))
529  num_degree = 15
530  end if
531 
532  ! find the triangle which (x,y) belongs to
533  ! i1, i2, i3 are the indices for the S1, S2 and S3 vertexes of the triangle
534  ! with S1 the vertex with the smallest y , S2 intermediate and S3 biggest
535  ! in other words , s1 bottom , s2 side, s3 top
536 
537  step = mesh%delta
538 
539  call sll_s_get_cell_vertices_index(x, y, mesh, i1, i2, i3)
540 
541  x1 = mesh%cartesian_coord(1, i1)
542  x2 = mesh%cartesian_coord(1, i2)
543  x3 = mesh%cartesian_coord(1, i3)
544  y1 = mesh%cartesian_coord(2, i1)
545  y2 = mesh%cartesian_coord(2, i2)
546  y3 = mesh%cartesian_coord(2, i3)
547 
548  k11 = mesh%hex_coord(1, i1)
549  k12 = mesh%hex_coord(2, i1)
550 
551  ! get the first 9 degrees of freedom
552 
553  ! values at the vertices of the triangle
554 
555  freedom(1) = f_tn(i1)
556  freedom(2) = f_tn(i2)
557  freedom(3) = f_tn(i3)
558 
559  ! values of the derivatives
560 
561  freedom(5) = deriv(3, i1) ! derivative from S1 to S3 (+h3)
562  freedom(8) = deriv(6, i3) ! derivative from S3 to S1 (-h3)
563 
564  test = .false.
565 
566  if (x2 <= x1) then ! first case : triangle oriented left
567 
568  freedom(4) = deriv(2, i1) ! derivative from S1 to S2 (+h2)
569  freedom(6) = deriv(5, i2) ! derivative from S2 to S1 (-h2)
570  freedom(7) = deriv(1, i2) ! derivative from S2 to S3 (+h1)
571  freedom(9) = deriv(4, i3) ! derivative from S3 to S2 (-h1)
572 
573  if (num_degree == 12) then
574  call get_normal_der(deriv, i1, i2, mesh, freedom(10:12))
575  elseif (num_degree == 15) then
576  call get_normal_der(deriv, i1, i2, mesh, freedom(13:15))
577  end if
578 
579  test = .true.
580 
581  else if (x2 > x1) then ! second case triangle oriented right
582  freedom(4) = deriv(1, i1) ! derivative from S1 to S2 (+h1)
583  freedom(6) = deriv(4, i2) ! derivative from S2 to S1 (-h1)
584  freedom(7) = deriv(2, i2) ! derivative from S2 to S3 (+h2)
585  freedom(9) = deriv(5, i3) ! derivative from S3 to S2 (-h2)
586 
587  if (num_degree == 12) then
588  call get_normal_der(deriv, i1, i2, mesh, freedom(10:12))
589  elseif (num_degree == 15) then
590  call get_normal_der(deriv, i1, i2, mesh, freedom(13:15))
591  end if
592 
593  test = .true.
594 
595  end if
596 
597  if (test .eqv. .false.) print *, "anomaly in sll_s_hermite_interpolation l289"
598 
599  a2 = 0.5_f64/aire
600  y3y = y3 - y
601  y2y = y2 - y
602  y1y = y1 - y
603  x3x = x3 - x
604  x2x = x2 - x
605  x1x = x1 - x
606 
607  l1 = a2*abs(x2x*y3y - x3x*y2y) ! barycentric coordinates
608  l2 = a2*abs(x1x*y3y - x3x*y1y)
609  l3 = 1._f64 - l1 - l2
610 
611  if (num_method == 9) then
612 
613  ! Computing the nine canonical basis functions
614  call base_zienkiewicz_9_degree_of_freedom(base, step, l1, l2, l3)
615 
616  else if (num_method == 10) then
617 
618  call sll_s_get_triangle_index(k11, k12, mesh, x, center_index)
619  freedom(10) = center_value(center_index)
620  !freedom(10) = (f_tn(i1)+f_tn(i2)+f_tn(i3))/3._f64
621  !call reconstruction_center_values(mesh,f_tn,center_index,freedom(10))
622 
623  ! Computing the ten canonical basis functions
624  call base_zienkiewicz_10_degree_of_freedom(base, step, l1, l2, l3)
625 
626  else if (num_method == 11) then
627 
628  ! Computing the basis for the cubic element of Hsieh-Clough-Tocher-reduced
629 
630  call base_hsieh_clough_tocher_reduced &
631  (base, x1, x3, y1, y2, y3, x, y, l1, l2, l3, mesh)
632 
633  else if (num_method == 12) then
634 
635  ! Computing the basis for the cubic element of Hsieh-Clough-Tocher-complete
636 
637  call base_hsieh_clough_tocher_complete &
638  (base, x1, x3, y1, y2, y3, x, y, l1, l2, l3, mesh)
639 
640  else if (num_method == 15) then
641 
642  ! values at the middle of the edges
643 
644  call sll_s_get_edge_index(k11, k12, mesh, x, edge_index1, edge_index2, edge_index3)
645  freedom(12) = edge_value(edge_index1)
646  freedom(11) = edge_value(edge_index2)
647  freedom(10) = edge_value(edge_index3)
648 
649  ! Computing the basis for the quartic element of Ganev_Dimitrov
650 
651  call base_ganev_dimitrov(base, l1, l2, l3, mesh)
652 
653  end if
654 
655  f = 0._f64
656 
657  do i = 1, num_degree
658  f = f + freedom(i)*base(i)
659  end do
660 
661  output_tn1(num) = f
662 
663  deallocate (freedom, base)
664 
665  end subroutine sll_s_hermite_interpolation
666 
667  subroutine base_zienkiewicz_9_degree_of_freedom(base, step, l1, l2, l3)
668  sll_real64, dimension(:), intent(out) :: base
669  sll_real64, intent(in) :: l1, l2, l3 ! barycentric coord
670  sll_real64, intent(in) :: step
671  sll_real64 :: phi, p05
672  sll_real64 :: l12, l22, l32, ksi1, ksi2, ksi3
673  sll_real64 :: ksi12, ksi13, ksi21, ksi23, ksi31, ksi32
674 
675  phi = l1*l2*l3 ! "bubble function"
676 
677  ksi1 = l1**3 - phi
678  ksi2 = l2**3 - phi
679  ksi3 = l3**3 - phi
680 
681  ! optimization variables
682 
683  p05 = phi*0.5_f64
684  l12 = l1**2
685  l22 = l2**2
686  l32 = l3**2
687 
688  ksi12 = l12*l2 + p05
689  ksi13 = l12*l3 + p05
690  ksi21 = l22*l1 + p05
691  ksi23 = l22*l3 + p05
692  ksi31 = l32*l1 + p05
693  ksi32 = l32*l2 + p05
694 
695  base(1) = 3._f64*l12 - 2._f64*ksi1
696  base(2) = 3._f64*l22 - 2._f64*ksi2
697  base(3) = 3._f64*l32 - 2._f64*ksi3
698  base(4) = step*ksi12
699  base(5) = step*ksi13
700  base(6) = step*ksi21
701  base(7) = step*ksi23
702  base(8) = step*ksi31
703  base(9) = step*ksi32
704 
705  end subroutine base_zienkiewicz_9_degree_of_freedom
706 
707  subroutine base_zienkiewicz_10_degree_of_freedom(base, step, l1, l2, l3)
708  sll_real64, dimension(:), intent(out) :: base
709  sll_real64, intent(in) :: l1, l2, l3 ! barycentric coord
710  sll_real64, intent(in) :: step
711  sll_real64 :: phi, p05, p9, p15
712  sll_real64 :: l12, l22, l32, ksi1, ksi2, ksi3
713  sll_real64 :: ksi12, ksi13, ksi21, ksi23, ksi31, ksi32
714 
715  phi = l1*l2*l3 ! "bubble function"
716 
717  ksi1 = l1**3 - phi
718  ksi2 = l2**3 - phi
719  ksi3 = l3**3 - phi
720 
721  ! optimization variables
722 
723  p05 = phi*0.5_f64
724  p9 = 9.0_f64*phi
725  p15 = 1.5_f64*phi
726  l12 = l1**2
727  l22 = l2**2
728  l32 = l3**2
729 
730  ksi12 = l12*l2 + p05
731  ksi13 = l12*l3 + p05
732  ksi21 = l22*l1 + p05
733  ksi23 = l22*l3 + p05
734  ksi31 = l32*l1 + p05
735  ksi32 = l32*l2 + p05
736 
737  base(1) = 3._f64*l12 - 2._f64*ksi1 - p9
738  base(2) = 3._f64*l22 - 2._f64*ksi2 - p9
739  base(3) = 3._f64*l32 - 2._f64*ksi3 - p9
740  base(4) = step*(ksi12 - p15)
741  base(5) = step*(ksi13 - p15)
742  base(6) = step*(ksi21 - p15)
743  base(7) = step*(ksi23 - p15)
744  base(8) = step*(ksi31 - p15)
745  base(9) = step*(ksi32 - p15)
746  base(10) = 27._f64*phi
747 
748  end subroutine base_zienkiewicz_10_degree_of_freedom
749 
750  subroutine get_sub_triangle_hex(mesh, x1, x3, y1, y2, y3, x, y, &
751  num_sub_triangle)
752  type(sll_t_hex_mesh_2d), pointer :: mesh
753  sll_real64, intent(in) :: x1, y1 ! cartesian coordinates of the lowest point
754  sll_real64, intent(in) :: y2 ! ordinate of the middle point
755  sll_real64, intent(in) :: x3, y3 ! cartesian coordinates of the highest point
756  sll_real64, intent(in) :: x, y ! cartesian coordinates of the point to be interpolated in a hexagonal mesh
757  sll_int32, intent(out) :: num_sub_triangle
758  sll_real64 :: slope
759 
760  slope = abs(mesh%r1_x1/mesh%r1_x2) ! sqrt(3._f64) -> default case
761 
762  ! what is the kind of the triangle T ?
763 
764  if (x < x1) then ! first kind : oriented left
765 
766  !finding the sub triangle K_l in which the interpolated point is
767 
768  if (y > y2) then! K1 or K2
769  if (y >= (x - x3)*slope + y3) then
770  num_sub_triangle = 1
771  else
772  num_sub_triangle = 2
773  end if
774  else ! K2 or K3
775  if (y >= -(x - x1)*slope + y1) then
776  num_sub_triangle = 2
777  else
778  num_sub_triangle = 3
779  end if
780  end if
781 
782  else ! second kind : oriented right
783 
784  !finding the triangle K_l in which the interpolated point is
785 
786  if (y > y2) then ! K1 or K2
787  if (y >= -(x - x3)*slope + y3) then
788  num_sub_triangle = 1
789  else
790  num_sub_triangle = 2
791  end if
792  else ! K2 or K3
793  if (y >= (x - x1)*slope + y1) then
794  num_sub_triangle = 2
795  else
796  num_sub_triangle = 3
797  end if
798  end if
799 
800  end if
801 
802  end subroutine get_sub_triangle_hex
803 
804  !*******************************************************************************
805  ! subroutines for hctr
806  !*******************************************************************************
807 
808  subroutine product_with_sigma_hct_r(li, lj, lk, xi)
809  sll_real64, dimension(:), intent(out) :: xi ! local base
810  sll_real64, intent(in) :: li, lj, lk ! barycentric coord
811  ! optimisation variables
812  sll_real64 :: li2j, li2k, lj2i
813  sll_real64 :: lj2k, lk2i, lk2j
814  sll_real64 :: li3, lj3, lk3
815  sll_real64 :: lijk
816 
817  li2k = li**2*lk
818  li2j = li**2*lj
819  lj2i = lj**2*li
820  lj2k = lj**2*lk
821  lk2i = lk**2*li
822  lk2j = lk**2*lj
823  li3 = li**3
824  lj3 = lj**3
825  lk3 = lk**3
826  lijk = li*lj*lk
827 
828  xi(1) = 4.5_f64*(li2k + li2j)
829  xi(2) = 0.5_f64*li3 + lj3 - 1.5_f64*li2k + 3.0_f64*(lj2i + lj2k + lijk)
830  xi(3) = 0.5_f64*li3 + lk3 - 1.5_f64*li2j + 3.0_f64*(lk2j + lk2i + lijk)
831  xi(4) = -0.25_f64*li3 + 1.25_f64*li2k + 0.5_f64*li2j
832  xi(5) = -0.25_f64*li3 + 0.5_f64*li2k + 1.25_f64*li2j
833  xi(6) = 0.25_f64*li3 - 0.5_f64*li2k - 0.25_f64*li2j + lj2i + lijk
834  xi(7) = -0.25_f64*li2k + 0.25_f64*li2j + lj2k + 0.5_f64*lijk
835  xi(8) = 0.25_f64*li2k - 0.25_f64*li2j + lk2j + 0.5_f64*lijk
836  xi(9) = 0.25_f64*li3 - 0.25_f64*li2k - 0.5_f64*li2j + lk2i + lijk
837 
838  end subroutine product_with_sigma_hct_r
839 
840  subroutine base_from_local_base_xi_htc_r(base, xi, num_sub_triangle, step)
841  sll_real64, dimension(:), intent(out) :: base
842  sll_real64, dimension(:), intent(in) :: xi
843  sll_real64, intent(in) :: step
844  sll_int32, intent(in) :: num_sub_triangle
845 
846  if (num_sub_triangle == 1) then
847  base(1) = xi(1)
848  base(2) = xi(2)
849  base(3) = xi(3)
850  base(4) = xi(5)*step
851  base(5) = xi(4)*step
852  base(6) = xi(6)*step
853  base(7) = xi(7)*step
854  base(8) = xi(9)*step
855  base(9) = xi(8)*step
856  elseif (num_sub_triangle == 2) then
857  base(1) = xi(3)
858  base(2) = xi(1)
859  base(3) = xi(2)
860  base(4) = xi(9)*step
861  base(5) = xi(8)*step
862  base(6) = xi(4)*step
863  base(7) = xi(5)*step
864  base(8) = xi(7)*step
865  base(9) = xi(6)*step
866  elseif (num_sub_triangle == 3) then
867  base(1) = xi(2)
868  base(2) = xi(3)
869  base(3) = xi(1)
870  base(4) = xi(7)*step
871  base(5) = xi(6)*step
872  base(6) = xi(8)*step
873  base(7) = xi(9)*step
874  base(8) = xi(5)*step
875  base(9) = xi(4)*step
876  end if
877 
878  end subroutine base_from_local_base_xi_htc_r
879 
880  subroutine base_hsieh_clough_tocher_reduced(base, x1, x3, y1, y2, y3, x, y, l1, l2, l3, mesh)
881  type(sll_t_hex_mesh_2d), pointer :: mesh
882  sll_real64, dimension(:), intent(out) :: base
883  sll_real64, intent(in) :: x1, x3, y1, y2, y3
884  sll_real64, intent(in) :: x, y, l1, l2, l3
885  sll_real64, dimension(:), allocatable :: xi
886  sll_real64 :: li, lj, lk, step
887  sll_int32 :: num_sub_triangle
888 
889  step = mesh%delta
890 
891  !finding the sub triangle K_l in which the interpolated point is
892 
893  call get_sub_triangle_hex(mesh, x1, x3, y1, y2, y3, x, y, &
894  num_sub_triangle)
895 
896  allocate (xi(1:9))
897 
898  if (num_sub_triangle == 1) then
899  li = l1
900  lj = l2
901  lk = l3
902  else if (num_sub_triangle == 2) then
903  li = l2
904  lj = l3
905  lk = l1
906  else if (num_sub_triangle == 3) then
907  li = l3
908  lj = l1
909  lk = l2
910  else
911  print *, "problem with finding the sub triangle"
912  end if
913 
914  call product_with_sigma_hct_r(li, lj, lk, xi)
915  call base_from_local_base_xi_htc_r(base, xi, num_sub_triangle, step)
916 
917  deallocate (xi)
918 
919  end subroutine base_hsieh_clough_tocher_reduced
920 
921  !*******************************************************************************
922  ! subroutines for hctc
923  !*******************************************************************************
924 
925  subroutine product_with_sigma_hct_c(li, lj, lk, xi)
926  sll_real64, dimension(:), intent(out) :: xi ! local base
927  sll_real64, intent(in) :: li, lj, lk ! barycentric coord
928  ! optimisation variables
929  sll_real64 :: li2j, li2k, lj2i, lj2k, lk2i, lk2j
930  sll_real64 :: li3, lj3, lk3
931  sll_real64 :: lijk
932 
933  li2k = li**2*lk
934  li2j = li**2*lj
935  lj2i = lj**2*li
936  lj2k = lj**2*lk
937  lk2i = lk**2*li
938  lk2j = lk**2*lj
939  li3 = li**3
940  lj3 = lj**3
941  lk3 = lk**3
942  lijk = li*lj*lk
943 
944  xi(1) = 4.5_f64*(li2k + li2j)
945  xi(2) = 0.5_f64*li3 + lj3 - 1.5_f64*li2k + 3.0_f64*(lj2i + lj2k + lijk)
946  xi(3) = 0.5_f64*li3 + lk3 - 1.5_f64*li2j + 3.0_f64*(lk2j + lk2i + lijk)
947  xi(4) = -1._f64/12._f64*li3 + 1.75_f64*li2k - 0.5_f64*li2j
948  xi(5) = -1._f64/12._f64*li3 - 0.5_f64*li2k + 1.75_f64*li2j
949  xi(6) = -7._f64/12._f64*li3 + 0.5_f64*li2k + 1.25_f64*li2j + lj2i - lijk
950  xi(7) = 2._f64/3._f64*li3 - 0.75_f64*li2k - 1.25_f64*li2j + lj2k + 1.5_f64*lijk
951  xi(8) = 2._f64/3._f64*li3 - 1.25_f64*li2k - 0.75_f64*li2j + lk2j + 1.5_f64*lijk
952  xi(9) = -7._f64/12._f64*li3 + 1.25_f64*li2k + 0.5_f64*li2j + lk2i - lijk
953  xi(10) = 4._f64/3._f64*li3 - 2._f64*li2k - 2._f64*li2j + 4._f64*lijk
954  xi(11) = -2._f64/3._f64*li3 + 2._f64*li2k
955  xi(12) = -2._f64/3._f64*li3 + 2._f64*li2j
956 
957  end subroutine product_with_sigma_hct_c
958 
959  subroutine base_from_local_base_xi_htc_c(base, xi, num_sub_triangle, mesh)
960  type(sll_t_hex_mesh_2d), pointer :: mesh
961  sll_real64, dimension(:), intent(out) :: base
962  sll_real64, dimension(:), intent(in) :: xi
963  sll_int32, intent(in) :: num_sub_triangle
964  sll_real64 :: step_sq, step
965 
966  step = mesh%delta
967  step_sq = abs(mesh%r1_x1)!step*abs(mesh%r1_x1)*real(mesh%num_cells,f64)/mesh%radius!sqrt(3._f64)*0.5_f64
968 
969  if (num_sub_triangle == 1) then
970  base(1) = xi(1)
971  base(2) = xi(2)
972  base(3) = xi(3)
973  base(4) = xi(5)*step
974  base(5) = xi(4)*step
975  base(6) = xi(6)*step
976  base(7) = xi(7)*step
977  base(8) = xi(9)*step
978  base(9) = xi(8)*step
979  base(10) = xi(10)*step_sq
980  base(11) = xi(11)*step_sq
981  base(12) = xi(12)*step_sq
982  elseif (num_sub_triangle == 2) then
983  base(1) = xi(3)
984  base(2) = xi(1)
985  base(3) = xi(2)
986  base(4) = xi(9)*step
987  base(5) = xi(8)*step
988  base(6) = xi(4)*step
989  base(7) = xi(5)*step
990  base(8) = xi(7)*step
991  base(9) = xi(6)*step
992  base(10) = xi(12)*step_sq
993  base(11) = xi(10)*step_sq
994  base(12) = xi(11)*step_sq
995  elseif (num_sub_triangle == 3) then
996  base(1) = xi(2)
997  base(2) = xi(3)
998  base(3) = xi(1)
999  base(4) = xi(7)*step
1000  base(5) = xi(6)*step
1001  base(6) = xi(8)*step
1002  base(7) = xi(9)*step
1003  base(8) = xi(5)*step
1004  base(9) = xi(4)*step
1005  base(10) = xi(11)*step_sq
1006  base(11) = xi(12)*step_sq
1007  base(12) = xi(10)*step_sq
1008  end if
1009 
1010  end subroutine base_from_local_base_xi_htc_c
1011 
1012  subroutine base_hsieh_clough_tocher_complete(base, x1, x3, y1, y2, y3, x, y, l1, l2, l3, mesh)
1013  type(sll_t_hex_mesh_2d), pointer :: mesh
1014  sll_real64, dimension(:), intent(out) :: base
1015  sll_real64, intent(in) :: x1, x3, y1, y2, y3
1016  sll_real64, intent(in) :: x, y, l1, l2, l3
1017  sll_real64, dimension(:), allocatable :: xi
1018  sll_real64 :: li, lj, lk
1019  sll_int32 :: num_sub_triangle
1020 
1021  !finding the sub triangle K_l in which the interpolated point is
1022 
1023  call get_sub_triangle_hex(mesh, x1, x3, y1, y2, y3, x, y, &
1024  num_sub_triangle)
1025 
1026  allocate (xi(1:12))
1027 
1028  if (num_sub_triangle == 1) then
1029  li = l1
1030  lj = l2
1031  lk = l3
1032  else if (num_sub_triangle == 2) then
1033  li = l2
1034  lj = l3
1035  lk = l1
1036  else if (num_sub_triangle == 3) then
1037  li = l3
1038  lj = l1
1039  lk = l2
1040  else
1041  print *, "problem with finding the sub triangle"
1042  end if
1043 
1044  call product_with_sigma_hct_c(li, lj, lk, xi)
1045  call base_from_local_base_xi_htc_c(base, xi, num_sub_triangle, mesh)
1046 
1047  deallocate (xi)
1048 
1049  end subroutine base_hsieh_clough_tocher_complete
1050 
1051  !*******************************************************************************
1052  ! computing the normal derivative
1053  !*******************************************************************************
1054 
1055  subroutine get_normal_der(deriv, i1, i2, mesh, freedom)
1056  type(sll_t_hex_mesh_2d), pointer :: mesh
1057  sll_real64, dimension(:, :), intent(in) :: deriv
1058  sll_real64, dimension(3), intent(out) :: freedom
1059  sll_real64 :: x1, x2
1060  sll_int32, intent(in) :: i1, i2
1061  sll_real64 :: fi_1, fi, fi1, fi2, fi_2, fi3
1062  sll_int32 :: ni_1, ni, ni1, ni2, ni_2, ni3
1063  sll_int32 :: h1, h2, num_cells, h16
1064  sll_int32 :: h11, h12, h13, h21, h22, h23, h14, h24
1065  sll_int32 :: h26, h2_4, h2_5, h25, h1_5, h1_4, h15
1066  sll_int32 :: h1_1, h1_3, h2_1, h2_3, h1_2, h2_2
1067  sll_real64, dimension(2) :: n1_l, n2_l, n3_l, n1_r, n2_r, n3_r
1068  sll_real64 :: det, v1, v2, v3, v4
1069 
1070  det = mesh%r1_x1*mesh%r2_x2 - mesh%r1_x2*mesh%r2_x1
1071 
1072  v1 = (+mesh%r1_x2*mesh%r2_x2 + mesh%r1_x1*mesh%r2_x1)/det
1073  v2 = (-mesh%r1_x2*mesh%r1_x2 - mesh%r1_x1*mesh%r1_x1)/det
1074  v3 = (+mesh%r2_x2*mesh%r2_x2 + mesh%r2_x1*mesh%r2_x1)/det
1075  v4 = (-mesh%r2_x2*mesh%r1_x2 - mesh%r2_x1*mesh%r1_x1)/det
1076 
1077  n1_l = (/v1, v2/)
1078  n2_l = (/-(v1 + v3), -(v2 + v4)/)
1079  n3_l = (/v3, v4/)
1080 
1081  ! n1_l = (/-1._f64/sqrt(3._f64) , -2._f64/sqrt(3._f64)/)
1082  ! n2_l = (/-1._f64/sqrt(3._f64) , +1._f64/sqrt(3._f64)/)
1083  ! n3_l = (/+2._f64/sqrt(3._f64) , +1._f64/sqrt(3._f64)/)
1084 
1085  n1_r = -n3_l
1086  n2_r = -n2_l
1087  n3_r = -n1_l
1088 
1089  num_cells = mesh%num_cells
1090 
1091  x1 = mesh%cartesian_coord(1, i1)
1092  x2 = mesh%cartesian_coord(1, i2)
1093 
1094  h1 = mesh%hex_coord(1, i1)
1095  h2 = mesh%hex_coord(2, i1)
1096 
1097  ! optimization variables
1098  h1_1 = h1 - 1
1099  h1_2 = h1 - 2
1100  h1_3 = h1 - 3
1101  h1_4 = h1 - 4
1102  h1_5 = h1 - 5
1103  h2_1 = h2 - 1
1104  h2_2 = h2 - 2
1105  h2_3 = h2 - 3
1106  h2_4 = h2 - 4
1107  h2_5 = h2 - 5
1108  h11 = h1 + 1
1109  h12 = h1 + 2
1110  h13 = h1 + 3
1111  h14 = h1 + 4
1112  h15 = h1 + 5
1113  h16 = h1 + 6
1114  h21 = h2 + 1
1115  h22 = h2 + 2
1116  h23 = h2 + 3
1117  h24 = h2 + 4
1118  h25 = h2 + 5
1119  h26 = h2 + 6
1120 
1121  ! let us determine the configuration : is the triangle oriented left or right ?
1122 
1123  if (x1 > x2) then ! oriented left
1124 
1125  ! the first normal derivative is oriented normal to r1 and m1 is in [S2;S3]
1126 
1127  if (test_in(h13, h26, num_cells)) then
1128  fi_2 = 0._f64
1129  else
1130  ni_2 = mesh%hex_to_global(h13, h26)
1131  fi_2 = deriv(1, ni_2)*n1_l(1) + deriv(2, ni_2)*n1_l(2)
1132  end if
1133 
1134  if (test_in(h12, h24, num_cells)) then
1135  fi_1 = 0._f64
1136  else
1137  ni_1 = mesh%hex_to_global(h12, h24)
1138  fi_1 = deriv(1, ni_1)*n1_l(1) + deriv(2, ni_1)*n1_l(2)
1139  end if
1140 
1141  if (test_in(h11, h22, num_cells)) then
1142  fi = 0._f64
1143  else
1144  ni = mesh%hex_to_global(h11, h22)
1145  fi = deriv(1, ni)*n1_l(1) + deriv(2, ni)*n1_l(2)
1146  end if
1147 
1148  if (test_in(h1, h2, num_cells)) then
1149  fi1 = 0._f64
1150  else
1151  ni1 = mesh%hex_to_global(h1, h2)
1152  fi1 = deriv(1, ni1)*n1_l(1) + deriv(2, ni1)*n1_l(2)
1153  end if
1154 
1155  if (test_in(h1_1, h2_2, num_cells)) then
1156  fi2 = 0._f64
1157  else
1158  ni2 = mesh%hex_to_global(h1_1, h2_2)
1159  fi2 = deriv(1, ni2)*n1_l(1) + deriv(2, ni2)*n1_l(2)
1160  end if
1161 
1162  if (test_in(h1_2, h2_4, num_cells)) then
1163  fi3 = 0._f64
1164  else
1165  ni3 = mesh%hex_to_global(h1_2, h2_4)
1166  fi3 = deriv(1, ni3)*n1_l(1) + deriv(2, ni3)*n1_l(2)
1167  end if
1168 
1169  !freedom(1) = ( -fi_1 + 9._f64*(fi + fi1) - fi2 )/16._f64
1170  freedom(1) = (3._f64*(fi_2 + fi3) - 25._f64*(fi_1 + fi2) + &
1171  150._f64*(fi + fi1))/256._f64
1172 
1173  ! the second normal derivative is oriented normal to r3 and m2 is in [S1;S3]
1174 
1175  if (test_in(h13, h2_2, num_cells)) then
1176  fi_2 = 0._f64
1177  else
1178  ni_2 = mesh%hex_to_global(h13, h2_2)
1179  fi_2 = deriv(1, ni_2)*n2_l(1) + deriv(2, ni_2)*n2_l(2)
1180  end if
1181 
1182  if (test_in(h12, h2_1, num_cells)) then
1183  fi_1 = 0._f64
1184  else
1185  ni_1 = mesh%hex_to_global(h12, h2_1)
1186  fi_1 = deriv(1, ni_1)*n2_l(1) + deriv(2, ni_1)*n2_l(2)
1187  end if
1188 
1189  if (test_in(h11, h2, num_cells)) then
1190  fi = 0._f64
1191  else
1192  ni = mesh%hex_to_global(h11, h2)
1193  fi = deriv(1, ni)*n2_l(1) + deriv(2, ni)*n2_l(2)
1194  end if
1195 
1196  if (test_in(h1, h2, num_cells)) then
1197  fi1 = 0._f64
1198  else
1199  ni1 = mesh%hex_to_global(h1, h21)
1200  fi1 = deriv(1, ni1)*n2_l(1) + deriv(2, ni1)*n2_l(2)
1201  end if
1202 
1203  if (test_in(h1_1, h22, num_cells)) then
1204  fi2 = 0._f64
1205  else
1206  ni2 = mesh%hex_to_global(h1_1, h22)
1207  fi2 = deriv(1, ni2)*n2_l(1) + deriv(2, ni2)*n2_l(2)
1208  end if
1209 
1210  if (test_in(h1_2, h23, num_cells)) then
1211  fi3 = 0._f64
1212  else
1213  ni3 = mesh%hex_to_global(h1_2, h23)
1214  fi3 = deriv(1, ni3)*n2_l(1) + deriv(2, ni3)*n2_l(2)
1215  end if
1216 
1217  !freedom(2) = ( -fi_1 + 9._f64*(fi + fi1) - fi2 )/16._f64
1218  freedom(2) = (3._f64*(fi_2 + fi3) - 25._f64*(fi_1 + fi2) + &
1219  150._f64*(fi + fi1))/256._f64
1220 
1221  ! the third normal derivative is oriented normal to r2 and m3 is in [S1;S2]
1222 
1223  if (test_in(h1_5, h2_2, num_cells)) then
1224  fi_2 = 0._f64
1225  else
1226  ni_2 = mesh%hex_to_global(h1_5, h2_2)
1227  fi_2 = deriv(1, ni_2)*n3_l(1) + deriv(2, ni_2)*n3_l(2)
1228  end if
1229 
1230  if (test_in(h1_3, h2_1, num_cells)) then
1231  fi_1 = 0._f64
1232  else
1233  ni_1 = mesh%hex_to_global(h1_3, h2_1)
1234  fi_1 = deriv(1, ni_1)*n3_l(1) + deriv(2, ni_1)*n3_l(2)
1235  end if
1236 
1237  if (test_in(h1_1, h2, num_cells)) then
1238  fi = 0._f64
1239  else
1240  ni = mesh%hex_to_global(h1_1, h2)
1241  fi = deriv(1, ni)*n3_l(1) + deriv(2, ni)*n3_l(2)
1242  end if
1243 
1244  if (test_in(h11, h21, num_cells)) then
1245  fi1 = 0._f64
1246  else
1247  ni1 = mesh%hex_to_global(h11, h21)
1248  fi1 = deriv(1, ni1)*n3_l(1) + deriv(2, ni1)*n3_l(2)
1249  end if
1250 
1251  if (test_in(h13, h22, num_cells)) then
1252  fi2 = 0._f64
1253  else
1254  ni2 = mesh%hex_to_global(h13, h22)
1255  fi2 = deriv(1, ni2)*n3_l(1) + deriv(2, ni2)*n3_l(2)
1256  end if
1257 
1258  if (test_in(h15, h23, num_cells)) then
1259  fi3 = 0._f64
1260  else
1261  ni3 = mesh%hex_to_global(h15, h23)
1262  fi3 = deriv(1, ni3)*n3_l(1) + deriv(2, ni3)*n3_l(2)
1263  end if
1264 
1265  !freedom(3) = ( -fi_1 + 9._f64*(fi + fi1) - fi2 )/16._f64
1266  freedom(3) = (3._f64*(fi_2 + fi3) - 25._f64*(fi_1 + fi2) + &
1267  150._f64*(fi + fi1))/256._f64
1268 
1269  else ! oriented right
1270 
1271  ! the first normal derivative is oriented normal to r2 and m1 is in [S2;S3]
1272 
1273  if (test_in(h16, h23, num_cells)) then
1274  fi_2 = 0._f64
1275  else
1276  ni_2 = mesh%hex_to_global(h16, h23)
1277  fi_2 = deriv(1, ni_2)*n1_r(1) + deriv(2, ni_2)*n1_r(2)
1278  end if
1279 
1280  if (test_in(h14, h22, num_cells)) then
1281  fi_1 = 0._f64
1282  else
1283  ni_1 = mesh%hex_to_global(h14, h22)
1284  fi_1 = deriv(1, ni_1)*n1_r(1) + deriv(2, ni_1)*n1_r(2)
1285  end if
1286 
1287  if (test_in(h12, h21, num_cells)) then
1288  fi = 0._f64
1289  else
1290  ni = mesh%hex_to_global(h12, h21)
1291  fi = deriv(1, ni)*n1_r(1) + deriv(2, ni)*n1_r(2)
1292  end if
1293 
1294  if (test_in(h1, h2, num_cells)) then
1295  fi1 = 0._f64
1296  else
1297  ni1 = mesh%hex_to_global(h1, h2)
1298  fi1 = deriv(1, ni1)*n1_r(1) + deriv(2, ni1)*n1_r(2)
1299  end if
1300 
1301  if (test_in(h1_2, h2_1, num_cells)) then
1302  fi2 = 0._f64
1303  else
1304  ni2 = mesh%hex_to_global(h1_2, h2_1)
1305  fi2 = deriv(1, ni2)*n1_r(1) + deriv(2, ni2)*n1_r(2)
1306  end if
1307 
1308  if (test_in(h1_4, h2_2, num_cells)) then
1309  fi3 = 0._f64
1310  else
1311  ni3 = mesh%hex_to_global(h1_4, h2_2)
1312  fi3 = deriv(1, ni3)*n1_r(1) + deriv(2, ni3)*n1_r(2)
1313  end if
1314 
1315  !freedom(1) = ( -fi_1 + 9._f64*(fi + fi1) - fi2 )/16._f64
1316  freedom(1) = (3._f64*(fi_2 + fi3) - 25._f64*(fi_1 + fi2) + &
1317  150._f64*(fi + fi1))/256._f64
1318 
1319  ! the second normal derivative is oriented normal to r3 and m2 is in [S1;S3]
1320 
1321  if (test_in(h1_2, h23, num_cells)) then
1322  fi_2 = 0._f64
1323  else
1324  ni_2 = mesh%hex_to_global(h1_2, h23)
1325  fi_2 = deriv(1, ni_2)*n2_r(1) + deriv(2, ni_2)*n2_r(2)
1326  end if
1327 
1328  if (test_in(h1_1, h22, num_cells)) then
1329  fi_1 = 0._f64
1330  else
1331  ni_1 = mesh%hex_to_global(h1_1, h22)
1332  fi_1 = deriv(1, ni_1)*n2_r(1) + deriv(2, ni_1)*n2_r(2)
1333  end if
1334 
1335  if (test_in(h1, h21, num_cells)) then
1336  fi = 0._f64
1337  else
1338  ni = mesh%hex_to_global(h1, h21)
1339  fi = deriv(1, ni)*n2_r(1) + deriv(2, ni)*n2_r(2)
1340  end if
1341 
1342  if (test_in(h11, h2, num_cells)) then
1343  fi1 = 0._f64
1344  else
1345  ni1 = mesh%hex_to_global(h11, h2)
1346  fi1 = deriv(1, ni1)*n2_r(1) + deriv(2, ni1)*n2_r(2)
1347  end if
1348 
1349  if (test_in(h12, h2_1, num_cells)) then
1350  fi2 = 0._f64
1351  else
1352  ni2 = mesh%hex_to_global(h12, h2_1)
1353  fi2 = deriv(1, ni2)*n2_r(1) + deriv(2, ni2)*n2_r(2)
1354  end if
1355 
1356  if (test_in(h13, h2_2, num_cells)) then
1357  fi3 = 0._f64
1358  else
1359  ni3 = mesh%hex_to_global(h13, h2_2)
1360  fi3 = deriv(1, ni3)*n2_r(1) + deriv(2, ni3)*n2_r(2)
1361  end if
1362 
1363  !freedom(2) = ( -fi_1 + 9._f64*(fi + fi1) - fi2 )/16._f64! interpolation
1364  freedom(2) = (3._f64*(fi_2 + fi3) - 25._f64*(fi_1 + fi2) + &
1365  150._f64*(fi + fi1))/256._f64
1366 
1367  ! the third normal derivative is oriented normal to r2 and m3 is in [S1;S2]
1368 
1369  if (test_in(h1_2, h2_5, num_cells)) then
1370  fi_2 = 0._f64
1371  else
1372  ni_2 = mesh%hex_to_global(h1_2, h2_5)
1373  fi_2 = deriv(1, ni_2)*n3_r(1) + deriv(2, ni_2)*n3_r(2)
1374  end if
1375 
1376  if (test_in(h1_1, h2_3, num_cells)) then
1377  fi_1 = 0._f64
1378  else
1379  ni_1 = mesh%hex_to_global(h1_1, h2_3)
1380  fi_1 = deriv(1, ni_1)*n3_r(1) + deriv(2, ni_1)*n3_r(2)
1381  end if
1382 
1383  if (test_in(h1, h2_1, num_cells)) then
1384  fi = 0._f64
1385  else
1386  ni = mesh%hex_to_global(h1, h2_1)
1387  fi = deriv(1, ni)*n3_r(1) + deriv(2, ni)*n3_r(2)
1388  end if
1389 
1390  if (test_in(h11, h21, num_cells)) then
1391  fi1 = 0._f64
1392  else
1393  ni1 = mesh%hex_to_global(h11, h21)
1394  fi1 = deriv(1, ni1)*n3_r(1) + deriv(2, ni1)*n3_r(2)
1395  end if
1396 
1397  if (test_in(h12, h23, num_cells)) then
1398  fi2 = 0._f64
1399  else
1400  ni2 = mesh%hex_to_global(h12, h23)
1401  fi2 = deriv(1, ni2)*n3_r(1) + deriv(2, ni2)*n3_r(2)
1402  end if
1403 
1404  if (test_in(h13, h25, num_cells)) then
1405  fi3 = 0._f64
1406  else
1407  ni3 = mesh%hex_to_global(h13, h25)
1408  fi3 = deriv(1, ni3)*n3_r(1) + deriv(2, ni3)*n3_r(2)
1409  end if
1410 
1411  !freedom(3) = ( -fi_1 + 9._f64*(fi + fi1) - fi2 )/16._f64! interpolation
1412  freedom(3) = (3._f64*(fi_2 + fi3) - 25._f64*(fi_1 + fi2) + &
1413  150._f64*(fi + fi1))/256._f64
1414 
1415  end if
1416 
1417  end subroutine get_normal_der
1418 
1419  function test_in(h1, h2, num_cells) result(bool)
1420  sll_int32 :: h1, h2, num_cells
1421  logical :: bool
1422 
1423  bool = abs(h1) > num_cells .or. abs(h2) > num_cells .or. &
1424  (h1)*(h2) < 0 .and. (abs(h1) + abs(h2) > num_cells)
1425 
1426  end function test_in
1427 
1428  subroutine product_with_sigma_ganev_dimitrov(li, lj, lk, xi)
1429  sll_real64, dimension(:), intent(out) :: xi ! local base
1430  sll_real64, intent(in) :: li, lj, lk ! barycentric coord
1431  ! optimisation variables
1432  sll_real64 :: li2, lj2, lk2
1433  sll_real64 :: li3, lj3, lk3
1434  sll_real64 :: li4, lj4, lk4
1435  sll_real64 :: li3j, li3k, lj3i, lj3k, lk3i, lk3j
1436  sll_real64 :: li2j2, li2k2, lj2k2
1437  sll_real64 :: li2jk, lij2k, lijk2
1438 
1439  !( with i = 1, j = 2 et k = 3 )
1440 
1441  li2 = li**2
1442  lj2 = lj**2
1443  lk2 = lk**2
1444  li2j2 = li2*lj2
1445  li2k2 = li2*lk2
1446  lj2k2 = lj2*lk2
1447  li3 = li**3
1448  lj3 = lj**3
1449  lk3 = lk**3
1450  li4 = li2**2
1451  lj4 = lj2**2
1452  lk4 = lk2**2
1453  li3j = li3*lj
1454  li3k = li3*lk
1455  lj3i = lj3*li
1456  lj3k = lj3*lk
1457  lk3i = lk3*li
1458  lk3j = lk3*lj
1459  li2jk = li2*lj*lk
1460  lij2k = li*lj2*lk
1461  lijk2 = li*lj*lk2
1462 
1463  xi(1) = li4 + 4._f64*(li3k + li3j) - 5._f64*(li2k2 + li2j2) - 4._f64*li2jk
1464  xi(2) = lj4 + 4._f64*(lj3i + lj3k) - 5._f64*(lj2k2 + li2j2) - 4._f64*lij2k
1465  xi(3) = lk4 + 4._f64*(lk3j + lk3i) - 5._f64*(lj2k2 + li2k2) - 4._f64*lijk2
1466  xi(4) = li3k - li2k2 + 0.5_f64*(-li2jk - lij2k + lijk2)
1467  xi(5) = li3j - li2j2 + 0.5_f64*(-li2jk + lij2k - lijk2)
1468  xi(6) = lj3i - li2j2 + 0.5_f64*(+li2jk - lij2k - lijk2)
1469  xi(7) = lj3k - lj2k2 + 0.5_f64*(-li2jk - lij2k + lijk2)
1470  xi(8) = lk3j - lj2k2 + 0.5_f64*(-li2jk + lij2k - lijk2)
1471  xi(9) = lk3i - li2k2 + 0.5_f64*(+li2jk - lij2k - lijk2)
1472  xi(10) = 16._f64*(lj2k2 - li2jk + lij2k + lijk2)
1473  xi(11) = 16._f64*(li2k2 + li2jk - lij2k + lijk2)
1474  xi(12) = 16._f64*(li2j2 + li2jk + lij2k - lijk2)
1475  xi(13) = 4._f64*(-li2jk + lij2k + lijk2)
1476  xi(14) = 4._f64*(li2jk - lij2k + lijk2)
1477  xi(15) = 4._f64*(li2jk + lij2k - lijk2)
1478 
1479  end subroutine product_with_sigma_ganev_dimitrov
1480 
1481  subroutine base_ganev_dimitrov(base, l1, l2, l3, mesh)
1482  type(sll_t_hex_mesh_2d), pointer :: mesh
1483  sll_real64, dimension(:), intent(out) :: base
1484  sll_real64, intent(in) :: l1, l2, l3
1485  sll_real64, dimension(:), allocatable :: xi
1486  sll_real64 :: step, step_sq
1487 
1488  step = mesh%delta
1489  step_sq = abs(mesh%r1_x1)!step*abs(mesh%r1_x1)*real(mesh%num_cells,f64)/mesh%radius!sqrt(3._f64)*0.5_f64
1490 
1491  allocate (xi(1:15))
1492 
1493  call product_with_sigma_ganev_dimitrov(l1, l2, l3, xi)
1494 
1495  base(1) = xi(1)
1496  base(2) = xi(2)
1497  base(3) = xi(3)
1498  base(4) = xi(5)*step
1499  base(5) = xi(4)*step
1500  base(6) = xi(6)*step
1501  base(7) = xi(7)*step
1502  base(8) = xi(9)*step
1503  base(9) = xi(8)*step
1504  base(10) = xi(10)
1505  base(11) = xi(11)
1506  base(12) = xi(12)
1507  base(13) = xi(13)*step_sq
1508  base(14) = xi(14)*step_sq
1509  base(15) = xi(15)*step_sq
1510 
1511  deallocate (xi)
1512 
1513  end subroutine base_ganev_dimitrov
1514 
1515  subroutine reconstruction_center_values(mesh, f_tn, center_index, center_value)
1516  type(sll_t_hex_mesh_2d), pointer :: mesh
1517  sll_real64, dimension(:), intent(in) :: f_tn
1518  sll_int32, intent(in):: center_index
1519  sll_real64, intent(out):: center_value
1520  sll_real64 :: center_value1, center_value2
1521  sll_real64 :: center_value3
1522  sll_int32 :: i1, i2, i3, k11, k12
1523  sll_int32 :: ni_3, ni_2, ni_1, ni, ni1, ni2, ni3, ni4
1524  sll_real64 :: x, y, x1
1525  !sll_real64 :: l1,l2,l3,l4,l5,l6,l7,l8
1526  sll_real64 :: r1, r2, r3, r4, r5, r6, r7, r8
1527  sll_real64 :: fi_3, fi_2, fi_1, fi, fi1, fi2, fi3, fi4
1528 
1529  x = mesh%center_cartesian_coord(1, center_index)
1530  y = mesh%center_cartesian_coord(2, center_index)
1531 
1532  call sll_s_get_cell_vertices_index(x, y, mesh, i1, i2, i3)
1533 
1534  x1 = mesh%cartesian_coord(1, i1)
1535 
1536  k11 = mesh%hex_coord(1, i1)
1537  k12 = mesh%hex_coord(2, i1)
1538 
1539  if (test_in(k11 - 3, k12 + 4, mesh%num_cells)) then
1540  fi_3 = 0._f64
1541  else
1542  ni_3 = mesh%hex_to_global(k11 - 3, k12 + 4)
1543  fi_3 = f_tn(ni_3)
1544  end if
1545 
1546  if (test_in(k11 - 2, k12 + 3, mesh%num_cells)) then
1547  fi_2 = 0._f64
1548  else
1549  ni_2 = mesh%hex_to_global(k11 - 2, k12 + 3)
1550  fi_2 = f_tn(ni_2)
1551  end if
1552 
1553  if (test_in(k11 - 1, k12 + 2, mesh%num_cells)) then
1554  fi_1 = 0._f64
1555  else
1556  ni_1 = mesh%hex_to_global(k11 - 1, k12 + 2)
1557  fi_1 = f_tn(ni_1)
1558  end if
1559 
1560  if (test_in(k11, k12 + 1, mesh%num_cells)) then
1561  fi = 0._f64
1562  else
1563  ni = mesh%hex_to_global(k11, k12 + 1)
1564  fi = f_tn(ni)
1565  end if
1566 
1567  if (test_in(k11 + 1, k12, mesh%num_cells)) then
1568  fi1 = 0._f64
1569  else
1570  ni1 = mesh%hex_to_global(k11 + 1, k12)
1571  fi1 = f_tn(ni1)
1572  end if
1573 
1574  if (test_in(k11 + 2, k12 - 1, mesh%num_cells)) then
1575  fi2 = 0._f64
1576  else
1577  ni2 = mesh%hex_to_global(k11 + 2, k12 - 1)
1578  fi2 = f_tn(ni2)
1579  end if
1580 
1581  if (test_in(k11 + 3, k12 - 2, mesh%num_cells)) then
1582  fi3 = 0._f64
1583  else
1584  ni3 = mesh%hex_to_global(k11 + 3, k12 - 2)
1585  fi3 = f_tn(ni3)
1586  end if
1587 
1588  if (test_in(k11 + 4, k12 - 3, mesh%num_cells)) then
1589  fi4 = 0._f64
1590  else
1591  ni4 = mesh%hex_to_global(k11 + 4, k12 - 3)
1592  fi4 = f_tn(ni4)
1593  end if
1594 
1595  ! l1 = 1.09739369e-2
1596  ! l2 = -9.60219478738e-2
1597  ! l3 = 0.768175583
1598  ! l4 = 0.3840877915
1599  ! l5 = -7.68175583e-2
1600  ! l6 = 9.60219478738e-3
1601 
1602  ! l1 = 0._f64
1603  ! l2 = -384._f64/6000._f64
1604  ! l3 = 1344._f64/2000._f64
1605  ! l4 = 896._f64/2000._f64
1606  ! l5 = -336._f64/6000._f64
1607  ! l6 = 0._f64
1608 
1609  ! l1 = -1.9704302961946860e-3_f64
1610  ! l2 = 1.9878164458669901e-2_f64
1611  ! l3 = -0.10671435656759629e0_f64
1612  ! l4 = 0.84482198949347154e0_f64
1613  ! l5 = 0.30720799617944411e0_f64
1614  ! l6 = -7.7983568260935790e-2_f64
1615  ! l7 = 1.6484331502311610e-2_f64
1616  ! l8 = -1.7241265091703488e-3_f64
1617 
1618  r8 = -1.9704302961946860e-3_f64
1619  r7 = 1.987816445866990e-2_f64
1620  r6 = -0.10671435656759629e0_f64
1621  r5 = 0.8448219894934715e0_f64
1622  r4 = 0.3072079961794441e0_f64
1623  r3 = -7.7983568260935790e-2_f64
1624  r2 = 1.6484331502311610e-2_f64
1625  r1 = -1.7241265091703488e-3_f64
1626 
1627  ! r1 = 0._f64
1628  ! r2 = 9.60219478738e-3_f64
1629  ! r3 = -7.68175583e-2_f64
1630  ! r4 = 0.3840877915_f64
1631  ! r5 = 0.768175583_f64
1632  ! r6 = -9.60219478738e-2_f64
1633  ! r7 = 1.09739369e-2_f64
1634  ! r8 = 0._f64
1635 
1636  ! r1 = 0._f64
1637  ! r2 = 0._f64
1638  ! r3 = -336._f64/6000._f64
1639  ! r4 = 896._f64/2000._f64
1640  ! r5 = 1344._f64/2000._f64
1641  ! r6 = -384._f64/6000._f64
1642  ! r7 = 0._f64
1643  ! r8 = 0._f64
1644 
1645  if (x < x1) then ! first case : triangle oriented left
1646  !center_value1 = l1*fi_3 + l2*fi_2 + l3*fi_1 + l4*fi + l5*fi1 + l6*fi2&
1647  ! + l7*fi3 + l8*fi4
1648  center_value1 = r8*fi_3 + r7*fi_2 + r6*fi_1 + r5*fi + r4*fi1 + r3*fi2 &
1649  + r2*fi3 + r1*fi4
1650  else ! second case : triangle oriented right
1651  center_value1 = r1*fi_3 + r2*fi_2 + r3*fi_1 + r4*fi + r5*fi1 + r6*fi2 &
1652  + r7*fi3 + r8*fi4
1653  end if
1654 
1655  if (x < x1) then ! first case : triangle oriented left
1656 
1657  if (test_in(k11 + 4, k12 + 8, mesh%num_cells)) then
1658  fi_3 = 0._f64
1659  else
1660  ni_3 = mesh%hex_to_global(k11 + 4, k12 + 8)
1661  fi_3 = f_tn(ni_3)
1662  end if
1663 
1664  if (test_in(k11 + 3, k12 + 6, mesh%num_cells)) then
1665  fi_2 = 0._f64
1666  else
1667  ni_2 = mesh%hex_to_global(k11 + 3, k12 + 6)
1668  fi_2 = f_tn(ni_2)
1669  end if
1670 
1671  if (test_in(k11 + 2, k12 + 4, mesh%num_cells)) then
1672  fi_1 = 0._f64
1673  else
1674  ni_1 = mesh%hex_to_global(k11 + 2, k12 + 4)
1675  fi_1 = f_tn(ni_1)
1676  end if
1677 
1678  if (test_in(k11 + 1, k12 + 2, mesh%num_cells)) then
1679  fi = 0._f64
1680  else
1681  ni = mesh%hex_to_global(k11 + 1, k12 + 2)
1682  fi = f_tn(ni)
1683  end if
1684 
1685  if (test_in(k11, k12, mesh%num_cells)) then
1686  fi1 = 0._f64
1687  else
1688  ni1 = mesh%hex_to_global(k11, k12)
1689  fi1 = f_tn(ni1)
1690  end if
1691 
1692  if (test_in(k11 - 1, k12 - 2, mesh%num_cells)) then
1693  fi2 = 0._f64
1694  else
1695  ni2 = mesh%hex_to_global(k11 - 1, k12 - 2)
1696  fi2 = f_tn(ni2)
1697  end if
1698 
1699  if (test_in(k11 - 2, k12 - 4, mesh%num_cells)) then
1700  fi3 = 0._f64
1701  else
1702  ni3 = mesh%hex_to_global(k11 - 2, k12 - 4)
1703  fi3 = f_tn(ni3)
1704  end if
1705 
1706  if (test_in(k11 - 3, k12 - 3, mesh%num_cells)) then
1707  fi4 = 0._f64
1708  else
1709  ni4 = mesh%hex_to_global(k11 - 3, k12 - 3)
1710  fi4 = f_tn(ni4)
1711  end if
1712 
1713  center_value2 = r1*fi_3 + r2*fi_2 + r3*fi_1 + r4*fi + r5*fi1 + r6*fi2 &
1714  + r7*fi3 + r8*fi4
1715 
1716  if (test_in(k11 - 7, k12 - 3, mesh%num_cells)) then
1717  fi_3 = 0._f64
1718  else
1719  ni_3 = mesh%hex_to_global(k11 - 7, k12 - 3)
1720  fi_3 = f_tn(ni_3)
1721  end if
1722 
1723  if (test_in(k11 - 5, k12 - 2, mesh%num_cells)) then
1724  fi_2 = 0._f64
1725  else
1726  ni_2 = mesh%hex_to_global(k11 - 5, k12 - 2)
1727  fi_2 = f_tn(ni_2)
1728  end if
1729 
1730  if (test_in(k11 - 3, k12 - 1, mesh%num_cells)) then
1731  fi_1 = 0._f64
1732  else
1733  ni_1 = mesh%hex_to_global(k11 - 3, k12 - 1)
1734  fi_1 = f_tn(ni_1)
1735  end if
1736 
1737  if (test_in(k11 - 1, k12, mesh%num_cells)) then
1738  fi = 0._f64
1739  else
1740  ni = mesh%hex_to_global(k11 - 1, k12)
1741  fi = f_tn(ni)
1742  end if
1743 
1744  if (test_in(k11 + 1, k12 + 1, mesh%num_cells)) then
1745  fi1 = 0._f64
1746  else
1747  ni1 = mesh%hex_to_global(k11 + 1, k12 + 1)
1748  fi1 = f_tn(ni1)
1749  end if
1750 
1751  if (test_in(k11 + 3, k12 + 2, mesh%num_cells)) then
1752  fi2 = 0._f64
1753  else
1754  ni2 = mesh%hex_to_global(k11 + 3, k12 + 2)
1755  fi2 = f_tn(ni2)
1756  end if
1757 
1758  if (test_in(k11 + 5, k12 + 3, mesh%num_cells)) then
1759  fi3 = 0._f64
1760  else
1761  ni3 = mesh%hex_to_global(k11 + 5, k12 + 3)
1762  fi3 = f_tn(ni3)
1763  end if
1764 
1765  if (test_in(k11 + 7, k12 + 4, mesh%num_cells)) then
1766  fi4 = 0._f64
1767  else
1768  ni4 = mesh%hex_to_global(k11 + 7, k12 + 4)
1769  fi4 = f_tn(ni4)
1770  end if
1771 
1772  center_value3 = r1*fi_3 + r2*fi_2 + r3*fi_1 + r4*fi + r5*fi1 + r6*fi2 &
1773  + r7*fi3 + r8*fi4
1774 
1775  else ! second case : triangle oriented right
1776 
1777  if (test_in(k11 + 4, k12 + 7, mesh%num_cells)) then
1778  fi_3 = 0._f64
1779  else
1780  ni_3 = mesh%hex_to_global(k11 + 4, k12 + 7)
1781  fi_3 = f_tn(ni_3)
1782  end if
1783 
1784  if (test_in(k11 + 3, k12 + 5, mesh%num_cells)) then
1785  fi_2 = 0._f64
1786  else
1787  ni_2 = mesh%hex_to_global(k11 + 3, k12 + 5)
1788  fi_2 = f_tn(ni_2)
1789  end if
1790 
1791  if (test_in(k11 + 2, k12 + 3, mesh%num_cells)) then
1792  fi_1 = 0._f64
1793  else
1794  ni_1 = mesh%hex_to_global(k11 + 2, k12 + 3)
1795  fi_1 = f_tn(ni_1)
1796  end if
1797 
1798  if (test_in(k11 + 1, k12 + 1, mesh%num_cells)) then
1799  fi = 0._f64
1800  else
1801  ni = mesh%hex_to_global(k11 + 1, k12 + 1)
1802  fi = f_tn(ni)
1803  end if
1804 
1805  if (test_in(k11, k12 - 1, mesh%num_cells)) then
1806  fi1 = 0._f64
1807  else
1808  ni1 = mesh%hex_to_global(k11, k12 - 1)
1809  fi1 = f_tn(ni1)
1810  end if
1811 
1812  if (test_in(k11 - 1, k12 - 3, mesh%num_cells)) then
1813  fi2 = 0._f64
1814  else
1815  ni2 = mesh%hex_to_global(k11 - 1, k12 - 3)
1816  fi2 = f_tn(ni2)
1817  end if
1818 
1819  if (test_in(k11 - 3, k12 - 5, mesh%num_cells)) then
1820  fi3 = 0._f64
1821  else
1822  ni3 = mesh%hex_to_global(k11 - 3, k12 - 5)
1823  fi3 = f_tn(ni3)
1824  end if
1825 
1826  if (test_in(k11 - 4, k12 - 7, mesh%num_cells)) then
1827  fi4 = 0._f64
1828  else
1829  ni4 = mesh%hex_to_global(k11 - 4, k12 - 7)
1830  fi4 = f_tn(ni4)
1831  end if
1832 
1833  center_value2 = r8*fi_3 + r7*fi_2 + r6*fi_1 + r5*fi + r4*fi1 + r3*fi2 &
1834  + r2*fi3 + r1*fi4
1835 
1836  if (test_in(k11 - 6, k12 - 3, mesh%num_cells)) then
1837  fi_3 = 0._f64
1838  else
1839  ni_3 = mesh%hex_to_global(k11 - 6, k12 - 3)
1840  fi_3 = f_tn(ni_3)
1841  end if
1842 
1843  if (test_in(k11 - 4, k12 - 2, mesh%num_cells)) then
1844  fi_2 = 0._f64
1845  else
1846  ni_2 = mesh%hex_to_global(k11 - 4, k12 - 2)
1847  fi_2 = f_tn(ni_2)
1848  end if
1849 
1850  if (test_in(k11 - 2, k12 - 1, mesh%num_cells)) then
1851  fi_1 = 0._f64
1852  else
1853  ni_1 = mesh%hex_to_global(k11 - 2, k12 - 1)
1854  fi_1 = f_tn(ni_1)
1855  end if
1856 
1857  if (test_in(k11, k12, mesh%num_cells)) then
1858  fi = 0._f64
1859  else
1860  ni = mesh%hex_to_global(k11, k12)
1861  fi = f_tn(ni)
1862  end if
1863 
1864  if (test_in(k11 + 2, k12 + 1, mesh%num_cells)) then
1865  fi1 = 0._f64
1866  else
1867  ni1 = mesh%hex_to_global(k11 + 2, k12 + 1)
1868  fi1 = f_tn(ni1)
1869  end if
1870 
1871  if (test_in(k11 + 4, k12 + 2, mesh%num_cells)) then
1872  fi2 = 0._f64
1873  else
1874  ni2 = mesh%hex_to_global(k11 + 4, k12 + 2)
1875  fi2 = f_tn(ni2)
1876  end if
1877 
1878  if (test_in(k11 + 6, k12 + 3, mesh%num_cells)) then
1879  fi3 = 0._f64
1880  else
1881  ni3 = mesh%hex_to_global(k11 + 6, k12 + 3)
1882  fi3 = f_tn(ni3)
1883  end if
1884 
1885  if (test_in(k11 + 8, k12 + 4, mesh%num_cells)) then
1886  fi4 = 0._f64
1887  else
1888  ni4 = mesh%hex_to_global(k11 + 8, k12 + 4)
1889  fi4 = f_tn(ni4)
1890  end if
1891 
1892  center_value3 = r8*fi_3 + r7*fi_2 + r6*fi_1 + r5*fi + r4*fi1 + r3*fi2 &
1893  + r2*fi3 + r1*fi4
1894  end if
1895 
1896  center_value = (center_value1 + center_value2 + center_value3)/3._f64
1897  !center_value = center_value3
1898 
1899  end subroutine reconstruction_center_values
1900 
1901  subroutine sll_s_print_method(num_method)
1902  sll_int32, intent(in) :: num_method
1903 
1904  if (num_method == 9) then
1905  print *, ""
1906  print *, "*********************************"
1907  print *, " Zienkiewicz_9_degree_of_freedom "
1908  print *, "*********************************"
1909  print *, ""
1910  else if (num_method == 10) then
1911  print *, ""
1912  print *, "*********************************"
1913  print *, " Zienkiewicz_10_degree_of_freedom"
1914  print *, "*********************************"
1915  print *, ""
1916  else if (num_method == 11) then
1917  print *, ""
1918  print *, "*********************************"
1919  print *, " Hsieh_Clough_Tocher_reduced "
1920  print *, "*********************************"
1921  print *, ""
1922  else if (num_method == 12) then
1923  print *, ""
1924  print *, "*********************************"
1925  print *, " Hsieh_Clough_Tocher_complete "
1926  print *, "*********************************"
1927  print *, ""
1928  else if (num_method == 15) then
1929  print *, ""
1930  print *, "*********************************"
1931  print *, " quartic element of Ganev_Dimitrov "
1932  print *, "*********************************"
1933  print *, ""
1934  else
1935  print *, "specify another number correspoonding to a existing implemented method 9, 10, 11, 12 or 15"
1936  end if
1937 
1938  end subroutine sll_s_print_method
1939 
1940 end module sll_m_interpolation_hex_hermite
sll_m_hermite_interpolation_2d
Definition: sll_m_hermite_interpolation_2d.F90:1
sll_m_hermite_interpolation_2d::sll_s_compute_w_hermite
subroutine, public sll_s_compute_w_hermite(w, r, s)
Definition: sll_m_hermite_interpolation_2d.F90:255
sll_m_hexagonal_meshes
Definition: sll_m_hexagonal_meshes.F90:23
sll_m_hexagonal_meshes::sll_s_get_triangle_index
subroutine, public sll_s_get_triangle_index(k1, k2, mesh, x, triangle_index)
Given a mesh point and the first cartesian coordinate of a point (that is not on the mesh) we return ...
Definition: sll_m_hexagonal_meshes.F90:1201
sll_m_hexagonal_meshes::sll_s_get_cell_vertices_index
subroutine, public sll_s_get_cell_vertices_index(x, y, mesh, s1, s2, s3)
Returns indices of the edges of a given cell.
Definition: sll_m_hexagonal_meshes.F90:1133
sll_m_hexagonal_meshes::sll_s_get_edge_index
subroutine, public sll_s_get_edge_index(k1, k2, mesh, x, edge_index1, edge_index2, edge_index3)
<BRIEF_DESCRIPTION>
Definition: sll_m_hexagonal_meshes.F90:1313
sll_m_interpolation_hex_hermite
Definition: sll_m_interpolation_hex_hermite.F90:1
sll_m_interpolation_hex_hermite::get_normal_der
subroutine get_normal_der(deriv, i1, i2, mesh, freedom)
Definition: sll_m_interpolation_hex_hermite.F90:1056
sll_m_interpolation_hex_hermite::base_from_local_base_xi_htc_c
subroutine base_from_local_base_xi_htc_c(base, xi, num_sub_triangle, mesh)
Definition: sll_m_interpolation_hex_hermite.F90:960
sll_m_interpolation_hex_hermite::get_sub_triangle_hex
subroutine get_sub_triangle_hex(mesh, x1, x3, y1, y2, y3, x, y, num_sub_triangle)
Definition: sll_m_interpolation_hex_hermite.F90:752
sll_m_interpolation_hex_hermite::sll_s_print_method
subroutine, public sll_s_print_method(num_method)
Definition: sll_m_interpolation_hex_hermite.F90:1902
sll_m_interpolation_hex_hermite::product_with_sigma_hct_c
subroutine product_with_sigma_hct_c(li, lj, lk, xi)
Definition: sll_m_interpolation_hex_hermite.F90:926
sll_m_interpolation_hex_hermite::product_with_sigma_ganev_dimitrov
subroutine product_with_sigma_ganev_dimitrov(li, lj, lk, xi)
Definition: sll_m_interpolation_hex_hermite.F90:1429
sll_m_interpolation_hex_hermite::base_zienkiewicz_10_degree_of_freedom
subroutine base_zienkiewicz_10_degree_of_freedom(base, step, l1, l2, l3)
Definition: sll_m_interpolation_hex_hermite.F90:708
sll_m_interpolation_hex_hermite::reconstruction_center_values
subroutine reconstruction_center_values(mesh, f_tn, center_index, center_value)
Definition: sll_m_interpolation_hex_hermite.F90:1516
sll_m_interpolation_hex_hermite::base_from_local_base_xi_htc_r
subroutine base_from_local_base_xi_htc_r(base, xi, num_sub_triangle, step)
Definition: sll_m_interpolation_hex_hermite.F90:841
sll_m_interpolation_hex_hermite::sll_s_hermite_interpolation
subroutine, public sll_s_hermite_interpolation(num, x, y, f_tn, center_value, edge_value, output_tn1, mesh, deriv, aire, num_method)
Definition: sll_m_interpolation_hex_hermite.F90:496
sll_m_interpolation_hex_hermite::base_ganev_dimitrov
subroutine base_ganev_dimitrov(base, l1, l2, l3, mesh)
Definition: sll_m_interpolation_hex_hermite.F90:1482
sll_m_interpolation_hex_hermite::product_with_sigma_hct_r
subroutine product_with_sigma_hct_r(li, lj, lk, xi)
Definition: sll_m_interpolation_hex_hermite.F90:809
sll_m_interpolation_hex_hermite::sll_s_der_finite_difference
subroutine, public sll_s_der_finite_difference(f_tn, p, step, mesh, deriv)
Definition: sll_m_interpolation_hex_hermite.F90:28
sll_m_interpolation_hex_hermite::base_hsieh_clough_tocher_complete
subroutine base_hsieh_clough_tocher_complete(base, x1, x3, y1, y2, y3, x, y, l1, l2, l3, mesh)
Definition: sll_m_interpolation_hex_hermite.F90:1013
sll_m_interpolation_hex_hermite::test_in
logical function test_in(h1, h2, num_cells)
Definition: sll_m_interpolation_hex_hermite.F90:1420
sll_m_interpolation_hex_hermite::base_zienkiewicz_9_degree_of_freedom
subroutine base_zienkiewicz_9_degree_of_freedom(base, step, l1, l2, l3)
Definition: sll_m_interpolation_hex_hermite.F90:668
sll_m_interpolation_hex_hermite::base_hsieh_clough_tocher_reduced
subroutine base_hsieh_clough_tocher_reduced(base, x1, x3, y1, y2, y3, x, y, l1, l2, l3, mesh)
Definition: sll_m_interpolation_hex_hermite.F90:881
sll_m_hexagonal_meshes::sll_t_hex_mesh_2d
2d hexagonal mesh
Definition: sll_m_hexagonal_meshes.F90:63
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1