Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_hexagonal_meshes.F90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------
2 ! SELALIB
3 !------------------------------------------------------------------------------
4 ! MODULE: sll_m_hexagonal_meshes
5 !
6 ! DESCRIPTION:
23 module sll_m_hexagonal_meshes
24 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
25 #include "sll_memory.h"
26 #include "sll_assert.h"
27 #include "sll_errors.h"
28 #include "sll_working_precision.h"
29 
30  use sll_m_constants, only: &
31  sll_p_sqrt3, &
32  sll_p_epsilon_0
33 
34  use sll_m_meshes_base, only: &
35  sll_c_mesh_2d_base
36 
37  use sll_m_tri_mesh_xmf, only: &
38  sll_s_write_tri_mesh_xmf
39 
40  implicit none
41 
42  public :: &
43  sll_f_cart_to_hex1, &
44  sll_f_cart_to_hex2, &
45  sll_f_cells_to_origin, &
46  sll_f_change_elements_notation, &
47  sll_o_delete, &
48  sll_s_delete_hex_mesh_2d, &
49  sll_s_display_hex_mesh_2d, &
50  sll_s_get_cell_vertices_index, &
51  sll_s_get_edge_index, &
52  sll_s_get_triangle_index, &
53  sll_f_local_to_global, &
54  sll_f_new_hex_mesh_2d, &
55  sll_s_hex_mesh_2d_init, &
56  sll_t_hex_mesh_2d, &
57  sll_s_write_caid_files
58 
59  private
60 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61 
63  type, extends(sll_c_mesh_2d_base) :: sll_t_hex_mesh_2d
64  ! A hexagonal mesh (composed by equilateral triangles)
65  ! is defined by three directional vectors (r1, r2, r3)
66  ! the number of cells, the radius, and the coordinates of the center
67  sll_int32 :: num_cells
68  sll_int32 :: num_pts_tot
69  sll_int32 :: num_triangles
70  sll_int32 :: num_edges
71  sll_real64 :: radius
72  sll_real64 :: center_x1
73  sll_real64 :: center_x2
74  sll_real64 :: delta
75  sll_real64 :: r1_x1
76  sll_real64 :: r1_x2
77  sll_real64 :: r2_x1
78  sll_real64 :: r2_x2
79  sll_real64 :: r3_x1
80  sll_real64 :: r3_x2
81 
82  ! Matrix containing mesh points coordinates in cartesian coordinates :
83  sll_real64, pointer, dimension(:, :) :: cartesian_coord
84  ! Matrix containing mesh points coordinates in hex coordinates (integers) :
85  sll_int32, pointer, dimension(:, :) :: hex_coord
86  ! Matrix containg global indices arranged from lower corner of hexagon
87  ! and following the r2, then r1 direction
88  sll_int32, pointer, dimension(:) :: global_indices
89 
90  ! matrix containing the cartesian coordinates of the triangles' centers
91  sll_real64, pointer, dimension(:, :) :: center_cartesian_coord
92  ! matrix containing the index of the respective center
93  ! of the 2 triangles at the top of most points
94  sll_int32, pointer, dimension(:, :) :: center_index
95 
96  ! The following 2 tables are not always needed, to avoid useless allocation
97  ! and initialization we define a flag to know if they are required
98  sll_int32 :: extra_tables
99  sll_real64, pointer, dimension(:, :) :: &
100  edge_center_cartesian_coord
101  ! matrix containing the index of the respective center of the 2 triangles
102  ! at the top of most points
103  sll_int32, pointer, dimension(:, :) :: edge_center_index
104 
105  contains
106  procedure, pass(mesh) :: eta1_node => eta1_node_hex
107  procedure, pass(mesh) :: eta2_node => eta2_node_hex
108  procedure, pass(mesh) :: eta1_cell_one_arg => eta1_cell_hex
109  procedure, pass(mesh) :: eta1_cell_two_arg => eta1_cell_hex_two_arg
110  procedure, pass(mesh) :: eta2_cell_one_arg => eta2_cell_hex
111  procedure, pass(mesh) :: eta2_cell_two_arg => eta2_cell_hex_two_arg
112  procedure, pass(mesh) :: index_hex_to_global
113  procedure, pass(mesh) :: hex_to_global
114  procedure, pass(mesh) :: global_to_hex1
115  procedure, pass(mesh) :: global_to_hex2
116  procedure, pass(mesh) :: global_to_x1
117  procedure, pass(mesh) :: global_to_x2
118  procedure, pass(mesh) :: sll_f_cart_to_hex1
119  procedure, pass(mesh) :: sll_f_cart_to_hex2
120  procedure, pass(mesh) :: global_to_local
121  procedure, pass(mesh) :: sll_f_local_to_global
122  procedure, pass(mesh) :: local_hex_to_global
123  procedure, pass(mesh) :: cell_type
124  procedure, pass(mesh) :: write_hex_mesh_2d
125  procedure, pass(mesh) :: write_hex_mesh_mtv
126  procedure, pass(mesh) :: get_neighbours
127  procedure, pass(mesh) :: hex_to_aligned_pt
128  procedure, pass(mesh) :: hex_to_aligned_elmt
129  procedure, pass(mesh) :: ref_to_hex_elmt
130  procedure, pass(mesh) :: ref_to_aligned_elmt
131  procedure, pass(mesh) :: display => sll_s_display_hex_mesh_2d
132  procedure, pass(mesh) :: write_field_hex_mesh_xmf
133  procedure, pass(mesh) :: delete => sll_s_delete_hex_mesh_2d
134  end type sll_t_hex_mesh_2d
135 
136  type hex_mesh_2d_ptr
137  type(sll_t_hex_mesh_2d), pointer :: hm
138  end type hex_mesh_2d_ptr
139 
140  interface sll_o_delete
145  module procedure sll_s_delete_hex_mesh_2d
146  end interface sll_o_delete
147 
148 contains
149 
150  ! Definition of a fonction to test if an argument is present
151  ! if it is it will asign it to the object at the slot,
152  ! else it will take a default value
153 #define TEST_PRESENCE_AND_ASSIGN_VAL( obj, arg, slot, default_val ) \
154  if (present(arg)) then; \
155  obj%slot = arg; \
156  else; \
157  obj%slot = default_val; \
158  end if
159 
160  !-------------------------------------------------------------------------
179  function sll_f_new_hex_mesh_2d( &
180  num_cells, &
181  center_x1, &
182  center_x2, &
183  r11, &
184  r12, &
185  r21, &
186  r22, &
187  r31, &
188  r32, &
189  radius, &
190  EXTRA_TABLES) result(mesh)
191 
192  type(sll_t_hex_mesh_2d), pointer :: mesh
193  sll_int32, intent(in) :: num_cells
194  sll_real64, optional, intent(in) :: radius
195  sll_real64, optional, intent(in) :: center_x1
196  sll_real64, optional, intent(in) :: center_x2
197  sll_real64, optional, intent(in) :: r11, r12
198  sll_real64, optional, intent(in) :: r21, r22
199  sll_real64, optional, intent(in) :: r31, r32
200  sll_int32, optional, intent(in) :: extra_tables
201  sll_int32 :: ierr
202 
203  sll_allocate(mesh, ierr)
204 
205  call sll_s_hex_mesh_2d_init( &
206  mesh, &
207  num_cells, &
208  radius, &
209  center_x1, &
210  center_x2, &
211  r11, &
212  r12, &
213  r21, &
214  r22, &
215  r31, &
216  r32, &
217  extra_tables)
218 
219  end function sll_f_new_hex_mesh_2d
220 
221  !-------------------------------------------------------------------------
237  subroutine sll_s_hex_mesh_2d_init( &
238  mesh, &
239  num_cells, &
240  radius, &
241  center_x1, &
242  center_x2, &
243  r1_x1, &
244  r1_x2, &
245  r2_x1, &
246  r2_x2, &
247  r3_x1, &
248  r3_x2, &
249  EXTRA_TABLES)
250 
251  type(sll_t_hex_mesh_2d) :: mesh
252  sll_int32, intent(in) :: num_cells
253  sll_real64, optional, intent(in) :: radius
254  sll_real64, optional, intent(in) :: center_x1
255  sll_real64, optional, intent(in) :: center_x2
256  sll_real64, optional, intent(in) :: r1_x1, r1_x2
257  sll_real64, optional, intent(in) :: r2_x1, r2_x2
258  sll_real64, optional, intent(in) :: r3_x1, r3_x2
259  sll_int32, optional, intent(in) :: extra_tables
260  sll_int32 :: ierr
261  sll_int32 :: i, j, global
262  sll_real64 :: position_x1
263  sll_real64 :: position_x2
264  sll_int32 :: k1, k2
265  sll_int32 :: index_tab
266  ! variables for optmizing computing time :
267  sll_int32 :: num_cells_plus1
268  sll_int32 :: num_cells_plus2
269 
270  ! By default the hexagonal mesh is centered at the (0,0) point
271  test_presence_and_assign_val(mesh, center_x1, center_x1, 0.0_f64)
272  test_presence_and_assign_val(mesh, center_x2, center_x2, 0.0_f64)
273 
274  ! By default the hexagonal mesh has a radius of 1.
275  test_presence_and_assign_val(mesh, radius, radius, 1.0_f64)
276 
277  ! By default the hexagonal mesh has for generator vectors :
278  ! r1 = (r11, r12) = ( sqrt(3)/2, 1/2 )
279  ! r2 = (r21, r22) = (-sqrt(3)/2, 1/2 )
280  ! r3 = (r31, r32) = ( 0, 1 )
281  test_presence_and_assign_val(mesh, r1_x1, r1_x1, sll_p_sqrt3*0.5_f64)
282  test_presence_and_assign_val(mesh, r1_x2, r1_x2, 0.5_f64)
283  test_presence_and_assign_val(mesh, r2_x1, r2_x1, -sll_p_sqrt3*0.5_f64)
284  test_presence_and_assign_val(mesh, r2_x2, r2_x2, 0.5_f64)
285  test_presence_and_assign_val(mesh, r3_x1, r3_x1, 0.0_f64)
286  test_presence_and_assign_val(mesh, r3_x2, r3_x2, 1.0_f64)
287 
288  mesh%num_cells = num_cells
289  mesh%delta = mesh%radius/real(num_cells, f64)
290 
291  ! The formula is = 6*sum(num_cells)+1 which simplifies to :
292  mesh%num_pts_tot = 3*mesh%num_cells*(mesh%num_cells + 1) + 1
293 
294  ! resizing :
295  mesh%r1_x1 = mesh%r1_x1*mesh%delta
296  mesh%r1_x2 = mesh%r1_x2*mesh%delta
297  mesh%r2_x1 = mesh%r2_x1*mesh%delta
298  mesh%r2_x2 = mesh%r2_x2*mesh%delta
299  mesh%r3_x1 = mesh%r3_x1*mesh%delta
300  mesh%r3_x2 = mesh%r3_x2*mesh%delta
301 
302  if (mesh%radius <= 0.0_f64) then
303  print *, 'ERROR, initialize_hex_mesh_2d(): ', &
304  'Problem to construct the mesh 2d '
305  print *, 'because radius <= 0.'
306  stop
307  end if
308  if (mesh%num_cells <= 0) then
309  print *, 'ERROR, initialize_hex_mesh_2d(): ', &
310  'Problem to construct the mesh 2d '
311  print *, 'because num_cells <= 0.'
312  stop
313  end if
314 
315  ! Allocation and initialization of coordinate matrices
316  ! and conectivity matrix
317  sll_allocate(mesh%cartesian_coord(2, mesh%num_pts_tot), ierr)
318  sll_allocate(mesh%hex_coord(2, mesh%num_pts_tot), ierr)
319  sll_allocate(mesh%global_indices(mesh%num_pts_tot), ierr)
320  mesh%cartesian_coord(:, :) = 0._f64
321  mesh%hex_coord(:, :) = 0
322  mesh%global_indices(:) = -1
323 
324  ! Initializing coordinates of first mesh point (ie. center of hexagon)
325  mesh%cartesian_coord(1, 1) = mesh%center_x1
326  mesh%cartesian_coord(2, 1) = mesh%center_x2
327 
328  !Useful variables
329  num_cells_plus1 = num_cells + 1
330  num_cells_plus2 = num_cells + 2
331 
332  ! ---------------------------------------------------------------------
333  ! BEGIN MATRICES INITIALIZATION ---------------------------------------
334 
335  ! Initializing coordinates of first mesh point (ie. center of hexagon)
336  global = 1
337  position_x1 = mesh%center_x1
338  position_x2 = mesh%center_x2 ! variable containing current position
339  mesh%cartesian_coord(1, global) = position_x1
340  mesh%cartesian_coord(2, global) = position_x2
341 
342  do i = 1, num_cells ! variable following r1
343  ! Incrementation on r1 direction as we are going to the next hexagon
344  position_x1 = position_x1 + mesh%r1_x1
345  position_x2 = position_x2 + mesh%r1_x2
346 
347  ! We follow each hexagon edge :
348  ! First edge
349  do j = 1, i ! following r2, the number of points on edge = i
350  global = global + 1
351 
352  mesh%cartesian_coord(1, global) = position_x1
353  mesh%cartesian_coord(2, global) = position_x2
354 
355  mesh%hex_coord(1, global) = i
356  mesh%hex_coord(2, global) = j - 1
357 
358  position_x1 = position_x1 + mesh%r2_x1
359  position_x2 = position_x2 + mesh%r2_x2
360  end do
361 
362  ! Second edge
363  do j = 1, i ! following -r1
364  global = global + 1
365 
366  mesh%cartesian_coord(1, global) = position_x1
367  mesh%cartesian_coord(2, global) = position_x2
368 
369  mesh%hex_coord(1, global) = i - j + 1
370  mesh%hex_coord(2, global) = i
371 
372  position_x1 = position_x1 - mesh%r1_x1
373  position_x2 = position_x2 - mesh%r1_x2
374  end do
375 
376  ! Third edge
377  do j = 1, i ! following -r3
378  global = global + 1
379 
380  mesh%cartesian_coord(1, global) = position_x1
381  mesh%cartesian_coord(2, global) = position_x2
382 
383  mesh%hex_coord(1, global) = -j + 1
384  mesh%hex_coord(2, global) = i - j + 1
385 
386  position_x1 = position_x1 - mesh%r3_x1
387  position_x2 = position_x2 - mesh%r3_x2
388  end do
389 
390  ! Fourth edge
391  do j = 1, i ! following -r2
392  global = global + 1
393 
394  mesh%cartesian_coord(1, global) = position_x1
395  mesh%cartesian_coord(2, global) = position_x2
396 
397  mesh%hex_coord(1, global) = -i
398  mesh%hex_coord(2, global) = -j + 1
399 
400  position_x1 = position_x1 - mesh%r2_x1
401  position_x2 = position_x2 - mesh%r2_x2
402  end do
403 
404  ! Fifth edge
405  do j = 1, i ! following r1
406  global = global + 1
407 
408  mesh%cartesian_coord(1, global) = position_x1
409  mesh%cartesian_coord(2, global) = position_x2
410 
411  mesh%hex_coord(1, global) = -i + j - 1
412  mesh%hex_coord(2, global) = -i
413 
414  position_x1 = position_x1 + mesh%r1_x1
415  position_x2 = position_x2 + mesh%r1_x2
416  end do
417 
418  ! Sixth edge
419  do j = 1, i ! following r3
420  global = global + 1
421 
422  mesh%cartesian_coord(1, global) = position_x1
423  mesh%cartesian_coord(2, global) = position_x2
424 
425  mesh%hex_coord(1, global) = j - 1
426  mesh%hex_coord(2, global) = -i + j - 1
427 
428  position_x1 = position_x1 + mesh%r3_x1
429  position_x2 = position_x2 + mesh%r3_x2
430  end do
431  end do
432 
433  ! Filling the global_indices table.
434  do global = 1, mesh%num_pts_tot
435 
436  k1 = mesh%hex_coord(1, global)
437  k2 = mesh%hex_coord(2, global)
438 
439  ! in order to switch from hexa coordinates to the global index
440  ! we need a routine that will compute a unique number index_tab
441  ! from the hexa coordinates. this unique number will index the array
442  ! global_indices which contains the numerotation
443 
444  call index_hex_to_global(mesh, k1, k2, index_tab)
445 
446  mesh%global_indices(index_tab) = global
447 
448  end do
449 
450  mesh%num_triangles = 6*num_cells*num_cells
451  sll_allocate(mesh%center_cartesian_coord(2, mesh%num_triangles), ierr)
452  sll_allocate(mesh%center_index(2, mesh%num_pts_tot), ierr)
453  ! we compute the cartesian coordinates and the indices
454  ! of the centers of the triangles of the mesh
455  call init_center_points_triangle(mesh)
456 
457  ! if needed we can compute the cartesian coordinates and the indices
458  ! of the edges' center of each triangle
459  mesh%num_edges = 3*num_cells*(3*num_cells + 1)
460 
461  test_presence_and_assign_val(mesh, extra_tables, extra_tables, 0)
462  if (mesh%EXTRA_TABLES .eq. 1) then
463  sll_allocate(mesh%edge_center_cartesian_coord(2, mesh%num_edges), ierr)
464  sll_allocate(mesh%edge_center_index(3, mesh%num_pts_tot), ierr)
465  call init_edge_center_triangle(mesh)
466  end if
467 
468  ! ----------------------------------------- END MATRICES INITIALIZATION
469  ! ---------------------------------------------------------------------
470 
471  end subroutine sll_s_hex_mesh_2d_init
472 
473 !---------------------------------------------------------------------------
478  subroutine init_center_points_triangle(mesh)
479  class(sll_t_hex_mesh_2d) :: mesh
480  sll_int32 :: center_index, global
481  sll_real64 :: k1, k2
482  sll_real64 :: x1, x2, x3
483  sll_real64 :: y1, y2, y3
484  sll_real64 :: xx, yy, r1x1
485  sll_real64 :: jacob, k1c, k2c
486  logical :: inside
487 
488  jacob = mesh%r1_x1*mesh%r2_x2 - mesh%r2_x1*mesh%r1_x2
489 
490  mesh%center_cartesian_coord(:, :) = 0._f64
491  mesh%center_index(:, :) = -1
492 
493  center_index = 0
494 
495  r1x1 = mesh%r1_x1*real(mesh%num_cells, f64)
496 
497  do global = 1, mesh%num_pts_tot
498  ! almost each point is the base of a lozenge , thus two triangle
499  ! from which we get two center points
500 
501  k1 = real(mesh%hex_coord(1, global), f64)
502  k2 = real(mesh%hex_coord(2, global), f64)
503 
504  ! center point in the left triangle
505 
506  x1 = k1*mesh%r1_x1 + k2*mesh%r2_x1
507  x2 = k1*mesh%r1_x1 + (k2 + 1)*mesh%r2_x1
508  x3 = (k1 + 1)*mesh%r1_x1 + (k2 + 1)*mesh%r2_x1
509  y1 = k1*mesh%r1_x2 + k2*mesh%r2_x2
510  y2 = k1*mesh%r1_x2 + (k2 + 1)*mesh%r2_x2
511  y3 = (k1 + 1)*mesh%r1_x2 + (k2 + 1)*mesh%r2_x2
512 
513  xx = x2 + ((x3 + x1)*0.5_f64 - x2)*2.0_f64/3.0_f64
514  yy = y2 + ((y3 + y1)*0.5_f64 - y2)*2.0_f64/3.0_f64
515 
516  !test to check if the triangle on the left is inside
517 
518  inside = .true.
519 
520  k1c = abs(mesh%r2_x2*xx - mesh%r2_x1*yy)/abs(jacob)
521  k2c = abs(mesh%r1_x1*yy - mesh%r1_x2*xx)/abs(jacob)
522 
523  if (ceiling(k1c) > mesh%num_cells) inside = .false.
524  if (ceiling(k2c) > mesh%num_cells) inside = .false.
525  if (abs(xx) > r1x1) inside = .false.
526 
527  if (inside) then
528  center_index = center_index + 1
529  mesh%center_cartesian_coord(1, center_index) = xx
530  mesh%center_cartesian_coord(2, center_index) = yy
531  mesh%center_index(1, global) = center_index
532  end if
533 
534  ! center point in the right triangle
535  x1 = k1*mesh%r1_x1 + k2*mesh%r2_x1
536  x2 = (k1 + 1)*mesh%r1_x1 + k2*mesh%r2_x1
537  x3 = (k1 + 1)*mesh%r1_x1 + (k2 + 1)*mesh%r2_x1
538  y1 = k1*mesh%r1_x2 + k2*mesh%r2_x2
539  y2 = (k1 + 1)*mesh%r1_x2 + k2*mesh%r2_x2
540  y3 = (k1 + 1)*mesh%r1_x2 + (k2 + 1)*mesh%r2_x2
541 
542  xx = x2 + ((x3 + x1)*0.5_f64 - x2)*2.0_f64/3.0_f64
543  yy = y2 + ((y3 + y1)*0.5_f64 - y2)*2.0_f64/3.0_f64
544 
545  ! test to check if the triangle on the left is inside
546 
547  inside = .true.
548 
549  k1c = abs(mesh%r2_x2*xx - mesh%r2_x1*yy)/abs(jacob)
550  k2c = abs(mesh%r1_x1*yy - mesh%r1_x2*xx)/abs(jacob)
551 
552  if (ceiling(k1c) > mesh%num_cells) inside = .false.
553  if (ceiling(k2c) > mesh%num_cells) inside = .false.
554  if (abs(xx) > r1x1) inside = .false.
555 
556  if (inside) then
557  center_index = center_index + 1
558  mesh%center_cartesian_coord(1, center_index) = xx
559  mesh%center_cartesian_coord(2, center_index) = yy
560  mesh%center_index(2, global) = center_index
561  end if
562 
563  end do
564 
565  end subroutine init_center_points_triangle
566 
567 !---------------------------------------------------------------------------
572  subroutine init_edge_center_triangle(mesh)
573  class(sll_t_hex_mesh_2d) :: mesh
574  sll_int32 :: edge_index
575  sll_int32 :: global, num_cells
576  sll_int32 :: k1, k2, k1c, k2c
577  sll_real64 :: x1, x2_l, x2_r, x3, y1, y2_l, y2_r, y3, xx, yy
578  logical :: inside
579 
580  num_cells = mesh%num_cells
581 
582  mesh%edge_center_cartesian_coord(:, :) = 0._f64
583  mesh%edge_center_index(:, :) = -1
584  edge_index = 0
585 
586  do global = 1, mesh%num_pts_tot
587  ! almost each point is the base of a lozenge , thus two triangle
588  ! from which we get two center points
589 
590  k1 = mesh%hex_coord(1, global)
591  k2 = mesh%hex_coord(2, global)
592 
593  x1 = real(k1, f64)*mesh%r1_x1 + real(k2, f64)*mesh%r2_x1
594  x2_l = real(k1, f64)*mesh%r1_x1 + real((k2 + 1), f64)*mesh%r2_x1
595  x2_r = real((k1 + 1), f64)*mesh%r1_x1 + real(k2, f64)*mesh%r2_x1
596  x3 = real((k1 + 1), f64)*mesh%r1_x1 + real((k2 + 1), f64)*mesh%r2_x1
597  y1 = real(k1, f64)*mesh%r1_x2 + real(k2, f64)*mesh%r2_x2
598  y2_l = real(k1, f64)*mesh%r1_x2 + real((k2 + 1), f64)*mesh%r2_x2
599  y2_r = real((k1 + 1), f64)*mesh%r1_x2 + real(k2, f64)*mesh%r2_x2
600  y3 = real((k1 + 1), f64)*mesh%r1_x2 + real((k2 + 1), f64)*mesh%r2_x2
601 
602  !test to check if the left edge is inside
603 
604  inside = .true.
605 
606  k1c = k1
607  k2c = k2 + 1
608 
609  if (abs(k1c) > num_cells .or. abs(k2c) > num_cells .or. &
610  (k1c)*(k2c) < 0 .and. (abs(k1c) + abs(k2c) > num_cells)) inside = .false.
611 
612  if (inside) then
613  ! center point of the left edge
614  xx = (x1 + x2_l)*0.5_f64
615  yy = (y1 + y2_l)*0.5_f64
616  edge_index = edge_index + 1
617  mesh%edge_center_cartesian_coord(1, edge_index) = xx
618  mesh%edge_center_cartesian_coord(2, edge_index) = yy
619  mesh%edge_center_index(1, global) = edge_index
620  end if
621 
622  !test to check if the middle edge is inside
623 
624  inside = .true.
625 
626  k1c = k1 + 1
627  k2c = k2 + 1
628 
629  if (abs(k1c) > num_cells .or. abs(k2c) > num_cells .or. &
630  (k1c)*(k2c) < 0 .and. (abs(k1c) + abs(k2c) > num_cells)) inside = .false.
631 
632  if (inside) then
633  ! center point of the middle edge
634  xx = (x1 + x3)*0.5_f64
635  yy = (y1 + y3)*0.5_f64
636  edge_index = edge_index + 1
637  mesh%edge_center_cartesian_coord(1, edge_index) = xx
638  mesh%edge_center_cartesian_coord(2, edge_index) = yy
639  mesh%edge_center_index(2, global) = edge_index
640  end if
641 
642  !test to check if the right edge is inside
643 
644  inside = .true.
645 
646  k1c = k1 + 1
647  k2c = k2
648 
649  if (abs(k1c) > num_cells .or. abs(k2c) > num_cells .or. &
650  (k1c)*(k2c) < 0 .and. (abs(k1c) + abs(k2c) > num_cells)) inside = .false.
651 
652  if (inside) then
653  ! center point of the right edge
654  xx = (x1 + x2_r)*0.5_f64
655  yy = (y1 + y2_r)*0.5_f64
656  edge_index = edge_index + 1
657  mesh%edge_center_cartesian_coord(1, edge_index) = xx
658  mesh%edge_center_cartesian_coord(2, edge_index) = yy
659  mesh%edge_center_index(3, global) = edge_index
660  end if
661 
662  end do
663 
664  end subroutine init_edge_center_triangle
665 
677  subroutine index_hex_to_global(mesh, k1, k2, index_tab)
678  class(sll_t_hex_mesh_2d) :: mesh
679  sll_int32, intent(in) :: k1, k2
680  sll_int32, intent(out) :: index_tab
681  sll_int32 :: k
682  sll_int32 :: nk1, nk2
683  sll_int32 :: n0
684  sll_int32 :: num_cells
685  sll_int32 :: num_cells_plus1
686 
687  num_cells = mesh%num_cells
688  num_cells_plus1 = num_cells + 1
689 
690  ! nk1 is the number of points before the edge corresponding to k1
691  ! nk2 is the number of points on the edge k1 :
692  ! if k1<=0 the points are in [-num_cells,k2]
693  ! if k1>0 ... [-num_cells+k1,k2]
694 
695  if (k1 .le. 0) then
696  k = num_cells + k1
697  !this value is always an integer, floor avoids the transformation
698  nk1 = num_cells*k + floor(real(k*(k + 1))*0.5)
699  nk2 = k2 + num_cells_plus1
700  else
701  ! n0 is the total number of points from (-num_cells,-num_cells) to
702  ! ( 0, numcells)
703  n0 = num_cells**2 + floor(real(num_cells*num_cells_plus1)*0.5)
704  nk1 = n0 + k1*(2*num_cells + 1) - floor(real(k1*(k1 - 1))*0.5)
705  nk2 = k2 + num_cells_plus1 - k1
706  end if
707 
708  index_tab = nk1 + nk2
709 
710  end subroutine index_hex_to_global
711 
718  function eta1_node_hex(mesh, i, j) result(res)
719  ! The coordinates (i, j) correspond to the (r1, r2) basis
720  ! This function returns the 1st coordinate on the cartesian system
721  class(sll_t_hex_mesh_2d), intent(in) :: mesh
722  sll_int32, intent(in) :: i
723  sll_int32, intent(in) :: j
724  sll_real64 :: res
725 
726  res = mesh%r1_x1*real(i, f64) + mesh%r2_x1*real(j, f64) + mesh%center_x1
727 
728  end function eta1_node_hex
729 
736  function eta2_node_hex(mesh, i, j) result(res)
737  ! The coordinates (k1, k2) correspond to the (r1, r2) basis
738  ! This function the 2nd coordinate on the cartesian system
739  class(sll_t_hex_mesh_2d), intent(in) :: mesh
740  sll_int32, intent(in) :: i
741  sll_int32, intent(in) :: j
742  sll_real64 :: res
743 
744  res = mesh%r1_x2*real(i, f64) + mesh%r2_x2*real(j, f64) + mesh%center_x2
745  end function eta2_node_hex
746 
752  function eta1_cell_hex(mesh, cell_num) result(res)
753  ! The index cell_num corresponds to the index of triangle
754  ! This function returns the 1st coordinate on the cartesian system
755  ! of the center of the triangle at num_ele
756  class(sll_t_hex_mesh_2d), intent(in) :: mesh
757  sll_int32, intent(in) :: cell_num
758  sll_real64 :: res
759 
760  sll_assert(0 < cell_num .and. cell_num <= mesh%num_triangles)
761  res = mesh%center_cartesian_coord(1, cell_num)
762  end function eta1_cell_hex
763 
769  function eta2_cell_hex(mesh, cell_num) result(res)
770  ! The index num_ele corresponds to the index of triangle
771  ! This function returns the 2nd coordinate on the cartesian system
772  ! of the center of the triangle at num_ele
773  class(sll_t_hex_mesh_2d), intent(in) :: mesh
774  sll_int32, intent(in) :: cell_num
775  sll_real64 :: res
776 
777  sll_assert(0 < cell_num .and. cell_num <= mesh%num_triangles)
778  res = mesh%center_cartesian_coord(2, cell_num)
779  end function eta2_cell_hex
780 
781  !---------------------------------------------------------------------------
786  function eta1_cell_hex_two_arg(mesh, i, j) result(res)
787  class(sll_t_hex_mesh_2d), intent(in) :: mesh
788  sll_int32, intent(in) :: i, j
789  sll_real64 :: res
790 
791  res = real(i + j + mesh%num_cells, f64)
792  print *, "Error in eta1_cell() :"
793  print *, " function only works with ONE parameter (num_cell)"
794  stop
795  end function eta1_cell_hex_two_arg
796 
797  !---------------------------------------------------------------------------
802  function eta2_cell_hex_two_arg(mesh, i, j) result(res)
803  class(sll_t_hex_mesh_2d), intent(in) :: mesh
804  sll_int32, intent(in) :: i, j
805  sll_real64 :: res
806 
807  res = real(i + j + mesh%num_cells, f64)
808  print *, "Error in eta2_cell() :"
809  print *, " function only works with ONE parameter (num_cell)"
810  stop
811  end function eta2_cell_hex_two_arg
812 
813  !---------------------------------------------------------------------------
821  function sll_f_cells_to_origin(k1, k2) result(val)
822  sll_int32, intent(in) :: k1
823  sll_int32, intent(in) :: k2
824  sll_int32 :: val
825 
826  ! We compute the number of cells from point to center
827  if (k1*k2 .gt. 0) then
828  val = max(abs(k1), abs(k2))
829  else
830  val = abs(k1) + abs(k2)
831  end if
832 
833  end function sll_f_cells_to_origin
834 
835  !---------------------------------------------------------------------------
844  subroutine cell_type(mesh, num_ele, val)
845  class(sll_t_hex_mesh_2d) :: mesh
846  sll_int32, intent(in) :: num_ele
847  sll_int32, intent(out) :: val
848  sll_int32 :: k1
849  sll_int32 :: k2
850  sll_int32 :: num_hex
851  sll_int32 :: lower_index
852  sll_real64 :: x1
853  sll_real64 :: y1
854  sll_real64 :: x2
855 
856  !Initialization:
857  val = -1
858 
859  ! Getting center coordinates
860  x1 = mesh%center_cartesian_coord(1, num_ele)
861  y1 = mesh%center_cartesian_coord(2, num_ele)
862  ! Getting hexagonal coordinates
863  k1 = sll_f_cart_to_hex1(mesh, x1, y1)
864  k2 = sll_f_cart_to_hex2(mesh, x1, y1)
865  !Getting number of cells to origin:
866  num_hex = sll_f_cells_to_origin(k1, k2)
867 
868  ! If the cell is not on the boundary, we can determine the
869  ! type of cell only by pairity of the cell index. If it's
870  ! even then is of type II, if it's odd then it's of type I.
871  if ((num_hex .lt. mesh%num_cells) .and. (mesh%num_cells .ne. 1)) then
872  if (modulo(num_ele, 2) .eq. 1) then
873  val = 1
874  else
875  val = 2
876  end if
877  elseif (mesh%num_cells .eq. 1) then
878  if ((num_ele .eq. 1) .or. (num_ele .eq. 4) .or. (num_ele .eq. 6)) then
879  val = 1
880  else
881  val = 2
882  end if
883  else
884  ! we just see if the center of the cell is to the
885  ! right or left of the lower point
886  lower_index = hex_to_global(mesh, k1, k2)
887  x2 = mesh%cartesian_coord(1, lower_index)
888  if (x2 < x1) then
889  val = 2
890  elseif (x2 > x1) then
891  val = 1
892  end if
893  end if
894  end subroutine cell_type
895 
896  !---------------------------------------------------------
906  function hex_to_global(mesh, k1, k2) result(val)
907  class(sll_t_hex_mesh_2d) :: mesh
908  sll_int32, intent(in) :: k1
909  sll_int32, intent(in) :: k2
910  sll_int32 :: distance
911  sll_int32 :: index_tab
912  sll_int32 :: val
913 
914  distance = sll_f_cells_to_origin(k1, k2)
915 
916  ! Test if we are in domain
917  if (distance .le. mesh%num_cells) then
918 
919  call index_hex_to_global(mesh, k1, k2, index_tab)
920  val = mesh%global_indices(index_tab)
921 
922  else
923  val = -1
924  end if
925 
926  end function hex_to_global
927 
928  !---------------------------------------------------------
934  function global_to_hex1(mesh, index) result(k1)
935  class(sll_t_hex_mesh_2d) :: mesh
936  sll_int32 :: index
937  sll_int32 :: k1
938 
939  k1 = mesh%hex_coord(1, index)
940  end function global_to_hex1
941 
942  !---------------------------------------------------------
948  function global_to_hex2(mesh, index) result(k2)
949  class(sll_t_hex_mesh_2d) :: mesh
950  sll_int32 :: index
951  sll_int32 :: k2
952 
953  k2 = mesh%hex_coord(2, index)
954  end function global_to_hex2
955 
956  !---------------------------------------------------------
962  function global_to_x1(mesh, index) result(x1)
963  class(sll_t_hex_mesh_2d) :: mesh
964  sll_int32 :: index
965  sll_real64 :: x1
966 
967  x1 = mesh%cartesian_coord(1, index)
968  end function global_to_x1
969 
970  !---------------------------------------------------------
976  function global_to_x2(mesh, index) result(x2)
977  class(sll_t_hex_mesh_2d) :: mesh
978  sll_int32 :: index
979  sll_real64 :: x2
980 
981  x2 = mesh%cartesian_coord(2, index)
982  end function global_to_x2
983 
984  !---------------------------------------------------------
991  function sll_f_cart_to_hex1(mesh, x1, x2) result(k1)
992  class(sll_t_hex_mesh_2d) :: mesh
993  sll_real64 :: x1
994  sll_real64 :: x2
995  sll_int32 :: k1
996  sll_real64 :: jacob
997 
998  jacob = mesh%r1_x1*mesh%r2_x2 - mesh%r2_x1*mesh%r1_x2
999  k1 = floor((mesh%r2_x2*x1 - mesh%r2_x1*x2)/jacob)
1000  end function sll_f_cart_to_hex1
1001 
1002  !---------------------------------------------------------
1009  function sll_f_cart_to_hex2(mesh, x1, x2) result(k2)
1010  class(sll_t_hex_mesh_2d) :: mesh
1011  sll_real64 :: x1
1012  sll_real64 :: x2
1013  sll_int32 :: k2
1014  sll_real64 :: jacob
1015 
1016  jacob = mesh%r1_x1*mesh%r2_x2 - mesh%r2_x1*mesh%r1_x2
1017  k2 = floor((mesh%r1_x1*x2 - mesh%r1_x2*x1)/jacob)
1018  end function sll_f_cart_to_hex2
1019 
1020  !---------------------------------------------------------
1029  function global_to_local(mesh, ref_index, global) result(local)
1030  class(sll_t_hex_mesh_2d) :: mesh
1031  sll_int32 :: ref_index
1032  sll_int32 :: global
1033  sll_int32 :: k1_ref, k2_ref
1034  sll_int32 :: k1_glob, k2_glob
1035  sll_int32 :: local
1036 
1037  if ((ref_index .le. mesh%num_pts_tot) .and. (ref_index .gt. 0) &
1038  .and. (global .le. mesh%num_pts_tot) .and. (global .gt. 0)) then
1039 
1040  k1_ref = mesh%hex_coord(1, ref_index)
1041  k2_ref = mesh%hex_coord(2, ref_index)
1042  k1_glob = mesh%hex_coord(1, global)
1043  k2_glob = mesh%hex_coord(2, global)
1044 
1045  local = mesh%hex_to_global(k1_ref - k1_glob, k2_ref - k2_glob)
1046  else
1047  ! Out of domain
1048  local = -1
1049  end if
1050 
1051  end function global_to_local
1052 
1053  !---------------------------------------------------------
1057  ! local index local_index in the ref_index system
1058  ! (see gloval_index(...) and global_to_local(...) for conventions)
1059  ! ie. sll_f_local_to_global(1, i) = i
1063  function sll_f_local_to_global(mesh, ref_index, local) result(global)
1064  class(sll_t_hex_mesh_2d) :: mesh
1065  sll_int32 :: ref_index, local
1066  sll_int32 :: k1_ref, k2_ref
1067  sll_int32 :: k1_loc, k2_loc
1068  sll_int32 :: global
1069 
1070  if ((ref_index .le. mesh%num_pts_tot) .and. (ref_index .gt. 0) &
1071  .and. (local .le. mesh%num_pts_tot) .and. (local .gt. 0)) then
1072  k1_ref = mesh%global_to_hex1(ref_index)
1073  k2_ref = mesh%global_to_hex2(ref_index)
1074  k1_loc = mesh%global_to_hex1(local)
1075  k2_loc = mesh%global_to_hex2(local)
1076 
1077  global = mesh%hex_to_global(k1_ref + k1_loc, k2_ref + k2_loc)
1078  else
1079  ! Out of domain
1080  global = -1
1081  end if
1082 
1083  end function sll_f_local_to_global
1084 
1085  !---------------------------------------------------------
1089  ! local index (k1_ref, k2_ref) in hex coordanites in the ref_index system
1090  ! (see gloval_index(...) and global_to_local(...) for conventions)
1091  ! ie. sll_f_local_to_global(1, i) = i
1096  function local_hex_to_global(mesh, k1_ref, k2_ref, local) result(global)
1097  class(sll_t_hex_mesh_2d) :: mesh
1098  sll_int32 :: k1_ref, k2_ref
1099  sll_int32 :: k1_loc, k2_loc
1100  sll_int32 :: local
1101  sll_int32 :: global
1102 
1103  k1_loc = mesh%global_to_hex1(local)
1104  k2_loc = mesh%global_to_hex2(local)
1105 
1106  if (sll_f_cells_to_origin(k1_ref + k1_loc, k2_ref + k2_loc) .lt. mesh%num_cells) then
1107  global = mesh%hex_to_global(k1_ref + k1_loc, k2_ref + k2_loc)
1108  else
1109  ! Out of domain
1110  global = -1
1111  end if
1112  end function local_hex_to_global
1113 
1114  !---------------------------------------------------------------------------
1121  ! mesh%center_cartesian_coord(2,cell_index), &
1132  subroutine sll_s_get_cell_vertices_index(x, y, mesh, s1, s2, s3)
1133  sll_real64, intent(in) :: x
1134  sll_real64, intent(in) :: y
1135  class(sll_t_hex_mesh_2d), intent(in) :: mesh ! Was pointer (YG - 05.10.2015)
1136  sll_int32, intent(out) :: s1
1137  sll_int32, intent(out) :: s2
1138  sll_int32, intent(out) :: s3
1139 
1140  sll_real64 :: xi, radius, step
1141  sll_int32 :: num_cells
1142  sll_int32 :: i, j
1143 
1144  num_cells = mesh%num_cells
1145  radius = mesh%radius
1146  step = mesh%delta
1147 
1148  ! converting (x,y) to hexagonal coordinates
1149  ! find the lowest point in the lozenge that contains (x,y)
1150  i = sll_f_cart_to_hex1(mesh, x, y)
1151  j = sll_f_cart_to_hex2(mesh, x, y)
1152 
1153  ! coordinates of the vertices of the lozenge :
1154  !(/i,j/),(/i,j+1/),(/i+1,j/), (/i+1,j+1/)
1155 
1156  ! coordinate of the abscisse that parts the lozenge
1157  ! in two equilateral triangle
1158 
1159  xi = (real(i, f64) - real(j, f64))*step*sll_p_sqrt3*0.5_f64
1160 
1161  ! testing which triangle (x,y) is in, which gives us its vertices'
1162  ! coordinates
1163 
1164  if (x > xi) then
1165  s1 = hex_to_global(mesh, i, j)
1166  s2 = hex_to_global(mesh, i + 1, j)
1167  s3 = hex_to_global(mesh, i + 1, j + 1)
1168  else if (x < xi) then
1169  s1 = hex_to_global(mesh, i, j)
1170  s2 = hex_to_global(mesh, i, j + 1)
1171  s3 = hex_to_global(mesh, i + 1, j + 1)
1172  else if (x == xi) then
1173  if (x < 0) then
1174  s1 = hex_to_global(mesh, i, j)
1175  s2 = hex_to_global(mesh, i + 1, j)
1176  s3 = hex_to_global(mesh, i + 1, j + 1)
1177  elseif (x >= 0) then
1178  s1 = hex_to_global(mesh, i, j)
1179  s2 = hex_to_global(mesh, i, j + 1)
1180  s3 = hex_to_global(mesh, i + 1, j + 1)
1181  end if
1182  end if
1183 
1184  end subroutine sll_s_get_cell_vertices_index
1185 
1186  !---------------------------------------------------------------------------
1200  subroutine sll_s_get_triangle_index(k1, k2, mesh, x, triangle_index)
1201  sll_int32, intent(in) :: k1
1202  sll_int32, intent(in) :: k2
1203  class(sll_t_hex_mesh_2d), intent(in) :: mesh ! Was pointer (YG - 05.10.2015)
1204  sll_real64, intent(in) :: x !cartessian_abscisse_other_vertice
1205  sll_int32, intent(out) :: triangle_index
1206 
1207  sll_int32 :: global
1208 
1209  triangle_index = -1
1210 
1211  ! almost every point is the lowest point of a lozenge , i.e. 2 triangles
1212  ! we get therefore 2 indices per points
1213  ! in order to have the correct one we test in which triangle we are
1214 
1215  global = hex_to_global(mesh, k1, k2)
1216 
1217  if ((global .gt. 0) .and. (global .le. mesh%num_pts_tot)) then
1218  if (x < mesh%cartesian_coord(1, global)) then
1219  triangle_index = mesh%center_index(1, global) !left triangle
1220  else
1221  triangle_index = mesh%center_index(2, global) !right triangle
1222  end if
1223  end if
1224 
1225  end subroutine sll_s_get_triangle_index
1226 
1227  !-------------------------------------------------------------------------
1235  subroutine get_neighbours(mesh, cell_index, nei_1, nei_2, nei_3)
1236  class(sll_t_hex_mesh_2d), intent(in) :: mesh
1237  sll_int32, intent(in) :: cell_index
1238  sll_int32, intent(out) :: nei_1
1239  sll_int32, intent(out) :: nei_2
1240  sll_int32, intent(out) :: nei_3
1241 
1242  sll_real64 :: xc
1243  sll_real64 :: yc
1244  sll_int32 :: s1
1245  sll_int32 :: s2
1246  sll_int32 :: s3
1247  sll_real64 :: x1, y1
1248  sll_real64 :: x2, y2
1249  sll_real64 :: x3, y3
1250  sll_real64 :: x
1251  sll_real64 :: y
1252  sll_int32 :: k1
1253  sll_int32 :: k2
1254  sll_real64 :: coef
1255 
1256  ! Getting the cell's center coordinates:
1257  xc = mesh%center_cartesian_coord(1, cell_index)
1258  yc = mesh%center_cartesian_coord(2, cell_index)
1259 
1260  ! Getting the cell's vertices indices:
1261  call sll_s_get_cell_vertices_index(xc, yc, mesh, s1, s2, s3)
1262 
1263  ! Getting the vertex coordinates:
1264  x1 = mesh%cartesian_coord(1, s1)
1265  y1 = mesh%cartesian_coord(2, s1)
1266  x2 = mesh%cartesian_coord(1, s2)
1267  y2 = mesh%cartesian_coord(2, s2)
1268  x3 = mesh%cartesian_coord(1, s3)
1269  y3 = mesh%cartesian_coord(2, s3)
1270 
1271  ! Getting the neighbours centers coordinates by symmetry to the cell's edges
1272  ! First center (symmetry with P1-P2) :
1273  coef = 2._f64*((xc - x1)*(x2 - x1) + &
1274  (yc - y1)*(y2 - y1))/((x2 - x1)**2 + &
1275  (y2 - y1)**2)
1276  x = coef*(x2 - x1) - (xc - x1) + x1
1277  y = coef*(y2 - y1) - (yc - y1) + y1
1278  ! Getting its index :
1279  k1 = sll_f_cart_to_hex1(mesh, x, y)
1280  k2 = sll_f_cart_to_hex2(mesh, x, y)
1281  call sll_s_get_triangle_index(k1, k2, mesh, x, nei_1)
1282 
1283  ! Second center (symmetry with P2-P3) :
1284  coef = 2._f64*((xc - x2)*(x3 - x2) + &
1285  (yc - y2)*(y3 - y2))/((x3 - x2)**2 + &
1286  (y3 - y2)**2)
1287  x = coef*(x3 - x2) - (xc - x2) + x2
1288  y = coef*(y3 - y2) - (yc - y2) + y2
1289  ! Getting its index :
1290  k1 = sll_f_cart_to_hex1(mesh, x, y)
1291  k2 = sll_f_cart_to_hex2(mesh, x, y)
1292  call sll_s_get_triangle_index(k1, k2, mesh, x, nei_2)
1293 
1294  ! Third center (symmetry with P1-P3) :
1295  coef = 2._f64*((xc - x1)*(x3 - x1) + &
1296  (yc - y1)*(y3 - y1))/((x3 - x1)**2 + &
1297  (y3 - y1)**2)
1298  x = coef*(x3 - x1) - (xc - x1) + x1
1299  y = coef*(y3 - y1) - (yc - y1) + y1
1300  ! Getting its index :
1301  k1 = sll_f_cart_to_hex1(mesh, x, y)
1302  k2 = sll_f_cart_to_hex2(mesh, x, y)
1303  call sll_s_get_triangle_index(k1, k2, mesh, x, nei_3)
1304 
1305  end subroutine get_neighbours
1306 
1307 !---------------------------------------------------------------------------
1312  subroutine sll_s_get_edge_index(k1, k2, mesh, x, edge_index1, edge_index2, edge_index3)
1313  type(sll_t_hex_mesh_2d), pointer :: mesh
1314  sll_real64, intent(in) :: x !cartessian_abscisse_other_vertice
1315  sll_int32, intent(in) :: k1, k2
1316  sll_int32, intent(out) :: edge_index1, edge_index2, edge_index3
1317  sll_int32 :: global, global2
1318 
1319  ! in short :
1320  ! returns the three indices of the edge of the triangle which contains
1321  ! a point of abscisse x and which lowest vertex is of hexa. coordinate
1322  ! k1, k2
1323 
1324  ! almost every point is the lowest point of a lozenge , i.e. 3 edges
1325  ! we get therefore 3 indices per points
1326  ! but we want the three indices of the triangle which contains a point
1327  ! of abscisse x.
1328  ! we therefore test in which kind of triangle we are
1329  ! ( oriented left or right ) then we get directly the indices of the first
1330  ! two edges then we deduce from the kind of triangle which point is to the
1331  ! left ( respectively to the right ) and get the index of the third edge
1332 
1333  global = hex_to_global(mesh, k1, k2)
1334 
1335  if (x < mesh%cartesian_coord(1, global)) then
1336  edge_index1 = mesh%edge_center_index(1, global)
1337  edge_index2 = mesh%edge_center_index(2, global)
1338  global2 = hex_to_global(mesh, k1, k2 + 1)
1339  edge_index3 = mesh%edge_center_index(3, global2)
1340  else
1341  edge_index1 = mesh%edge_center_index(3, global)
1342  edge_index2 = mesh%edge_center_index(2, global)
1343  global2 = hex_to_global(mesh, k1 + 1, k2)
1344  edge_index3 = mesh%edge_center_index(1, global2)
1345  end if
1346 
1347  if (edge_index1 == -1) print *, "problem in sll_s_get_edge_index l", __line__
1348  if (edge_index2 == -1) print *, "problem in sll_s_get_edge_index l", __line__
1349  if (edge_index3 == -1) print *, "problem in sll_s_get_edge_index l", __line__
1350 
1351  end subroutine sll_s_get_edge_index
1352 
1370  function sll_f_change_elements_notation(mesh, i_elmt_old) result(i_elmt)
1371  type(sll_t_hex_mesh_2d), pointer :: mesh
1372  sll_int32, intent(in) :: i_elmt_old
1373  sll_int32 :: i_elmt
1374  sll_int32 :: e1, e2, e3
1375  sll_int32 :: j1, j2, j3
1376  sll_int32 :: e1_k1, e1_k2
1377  sll_int32 :: e2_k1, e2_k2
1378  sll_int32 :: e3_k1, e3_k2
1379  sll_int32 :: e1_distance
1380  sll_int32 :: e2_distance
1381  sll_int32 :: e3_distance
1382  sll_int32 :: sixth
1383  sll_int32 :: layer
1384  sll_int32 :: npts_layer
1385  sll_int32 :: npts_layer_1
1386  sll_int32 :: displacement
1387 
1388  if (i_elmt_old .eq. -1) then
1389  i_elmt = -1
1390  return
1391  end if
1392 
1393  ! Getting cell vertices
1394  call sll_s_get_cell_vertices_index(mesh%center_cartesian_coord(1, i_elmt_old), &
1395  mesh%center_cartesian_coord(2, i_elmt_old), &
1396  mesh, &
1397  e1, e2, e3)
1398 
1399  ! Getting their hexagonal coordinates
1400  e1_k1 = mesh%hex_coord(1, e1)
1401  e1_k2 = mesh%hex_coord(2, e1)
1402  e2_k1 = mesh%hex_coord(1, e2)
1403  e2_k2 = mesh%hex_coord(2, e2)
1404  e3_k1 = mesh%hex_coord(1, e3)
1405  e3_k2 = mesh%hex_coord(2, e3)
1406 
1407  ! Getting their distance to origin:
1408  e1_distance = sll_f_cells_to_origin(e1_k1, e1_k2)
1409  e2_distance = sll_f_cells_to_origin(e2_k1, e2_k2)
1410  e3_distance = sll_f_cells_to_origin(e3_k1, e3_k2)
1411 
1412  ! computing on which hexagonal layer is the cell
1413  layer = e1_distance + e2_distance + e3_distance
1414  layer = layer/3
1415 
1416  !Computing on which sixth of the hexagon are we:
1417  if ((e1_k1 >= 0) .and. (e1_k2 >= 0) .and. &
1418  (e2_k1 >= 0) .and. (e2_k2 >= 0) .and. &
1419  (e3_k1 >= 0) .and. (e3_k2 >= 0)) then
1420  if (e1_k1 + e2_k1 + e3_k1 .gt. e1_k2 + e2_k2 + e3_k2) then
1421  sixth = 1
1422  else
1423  sixth = 2
1424  end if
1425  elseif ((e1_k1 <= 0) .and. (e1_k2 <= 0) .and. &
1426  (e2_k1 <= 0) .and. (e2_k2 <= 0) .and. &
1427  (e3_k1 <= 0) .and. (e3_k2 <= 0)) then
1428  if (abs(e1_k1 + e2_k1 + e3_k1) .gt. abs(e1_k2 + e2_k2 + e3_k2)) then
1429  sixth = 4
1430  else
1431  sixth = 5
1432  end if
1433  elseif ((e1_k1 <= 0) .and. (e1_k2 >= 0) .and. &
1434  (e2_k1 <= 0) .and. (e2_k2 >= 0) .and. &
1435  (e3_k1 <= 0) .and. (e3_k2 >= 0)) then
1436  sixth = 3
1437  else
1438  sixth = 6
1439  end if
1440 
1441  if (layer .eq. 0) then
1442  ! Treating the first layer separetly
1443  i_elmt = sixth
1444  else
1445  ! founding the displacement:
1446  npts_layer = 3*(layer + 1)*layer + 1
1447  npts_layer_1 = 3*(layer - 1)*layer + 1
1448  j1 = e1 - npts_layer
1449  j2 = e2 - npts_layer
1450  j3 = e3 - npts_layer
1451  if (j1 <= 0) then
1452  if (j2 <= 0) then
1453  displacement = e3 - npts_layer
1454  if (mod(sixth, 2) .eq. 1) then
1455  displacement = 2*(displacement - 1) - (sixth - 1)
1456  else
1457  displacement = 2*(displacement - 1) - (sixth - 1)
1458  end if
1459  elseif (j3 <= 0) then
1460  displacement = e2 - npts_layer
1461  if (mod(sixth, 2) .eq. 1) then
1462  displacement = 2*(displacement - 1) - (sixth - 1)
1463  else
1464  displacement = 2*(displacement - 1) - (sixth - 1)
1465  end if
1466  else
1467  displacement = e1 - npts_layer_1
1468  if (mod(sixth, 2) .eq. 1) then
1469  displacement = 2*(displacement - 1) + sixth
1470  else
1471  displacement = 2*(displacement - 1) + sixth
1472  end if
1473  end if
1474  elseif (j2 <= 0) then
1475  if (j3 <= 0) then
1476  displacement = e1 - npts_layer
1477  if (mod(sixth, 2) .eq. 1) then
1478  displacement = 2*(displacement - 1) - (sixth - 1)
1479  else
1480  displacement = 2*(displacement - 1) - (sixth - 1)
1481  end if
1482  else
1483  displacement = e2 - npts_layer_1
1484  if ((sixth .eq. 6) .and. (j1 .eq. 1)) then
1485  displacement = e3 - npts_layer_1
1486  elseif ((sixth .eq. 6) .and. (j3 .eq. 1)) then
1487  displacement = e1 - npts_layer_1
1488  else
1489  displacement = 2*(displacement - 1) + sixth
1490  end if
1491  end if
1492  else
1493  displacement = e3 - npts_layer_1
1494  if (mod(sixth, 2) .eq. 1) then
1495  displacement = 2*(displacement - 1) + sixth
1496  else
1497  displacement = 2*(displacement - 1) + sixth
1498  end if
1499  end if
1500 
1501  ! result
1502  i_elmt = 6*layer*layer + displacement
1503  end if
1504  end function sll_f_change_elements_notation
1505 
1506  !---------------------------------------------------------------------------
1516  recursive subroutine hex_to_aligned_pt(mesh, ind, transf, x_new, y_new)
1517  class(sll_t_hex_mesh_2d), intent(in) :: mesh
1518  sll_int32, intent(in) :: ind
1519  character(len=*), intent(in) :: transf
1520  sll_real64, intent(out) :: x_new
1521  sll_real64, intent(out) :: y_new
1522  ! Local
1523  sll_real64 :: x
1524  sll_real64 :: y
1525  sll_real64 :: radius
1526  sll_real64 :: angle
1527  sll_real64 :: asin_gamma
1528  sll_real64 :: kappa
1529  sll_real64 :: x_temp
1530  sll_real64 :: y_temp
1531  sll_real64 :: dist_to_origin
1532  sll_int32 :: k1
1533  sll_int32 :: k2
1534  sll_int32 :: cells_dist
1535 
1536  ! Initalization
1537  x = mesh%cartesian_coord(1, ind)
1538  y = mesh%cartesian_coord(2, ind)
1539  x_new = -1000.0_f64
1540  y_new = -1000.0_f64
1541 
1542  select case (transf)
1543  case ("CIRCLE")
1544  !........................................................................
1545  ! Computing current radius from point to origin
1546  dist_to_origin = sqrt(x**2 + y**2)
1547 
1548  ! Computing the actual radius if the point was on a circle:
1549  ! First we get the hexagonal coordinates
1550  k1 = mesh%hex_coord(1, ind)
1551  k2 = mesh%hex_coord(2, ind)
1552  cells_dist = sll_f_cells_to_origin(k1, k2)
1553  radius = mesh%radius/real(mesh%num_cells, f64)*real(cells_dist, f64)
1554 
1555  if (dist_to_origin .le. sll_p_epsilon_0) then
1556  x_new = 0._f64
1557  y_new = 0._f64
1558  else
1559  x_new = x*radius/dist_to_origin
1560  y_new = y*radius/dist_to_origin
1561  end if
1562  !........................................................................
1563  case ("TOKAMAK")
1564  !........................................................................
1565  call mesh%hex_to_aligned_pt(ind, "CIRCLE", x_temp, y_temp)
1566  ! Computing current radius and angle
1567  radius = sqrt(x_temp**2 + y_temp**2)
1568  angle = 0._f64
1569  if (.not. ((x_temp .eq. 0._f64) .and. (y_temp .eq. 0._f64))) then
1570  angle = atan2(y_temp, x_temp)
1571  end if
1572 
1573  ! Some constants taken from "Noncircular, finite aspect ratio, local
1574  ! equilibrium model" by R.L. Miller et al, Physics of Plamas, Vol 5 ('98)
1575  asin_gamma = 0.4290421956533866_f64
1576  kappa = 1.66_f64
1577  x_new = radius*cos(angle + asin_gamma*sin(angle))
1578  y_new = kappa*radius*sin(angle)
1579  !........................................................................
1580  case default
1581  print *, "ERROR in hex_to_aligned_pt: no known transformation =>", transf
1582  print *, " Options are : CIRCLE and TOKAMAK."
1583  stop
1584  end select
1585  end subroutine hex_to_aligned_pt
1586 
1587  !---------------------------------------------------------------------------
1596  subroutine hex_to_aligned_elmt(mesh, i_elmt, transf, transf_matA, transf_vecB)
1597  class(sll_t_hex_mesh_2d), intent(in) :: mesh
1598  sll_int32, intent(in) :: i_elmt
1599  character(len=*), intent(in) :: transf
1600  sll_real64, dimension(2, 2), intent(out) :: transf_mata
1601  sll_real64, dimension(2), intent(out) :: transf_vecb
1602  ! Local
1603  sll_real64 :: x_v1, y_v1
1604  sll_real64 :: x_v2, y_v2
1605  sll_real64 :: x_v3, y_v3
1606  sll_real64 :: xp_v1, yp_v1
1607  sll_real64 :: xp_v2, yp_v2
1608  sll_real64 :: xp_v3, yp_v3
1609  sll_real64 :: ccx, ccy
1610  sll_int32 :: e1
1611  sll_int32 :: e2
1612  sll_int32 :: e3
1613  sll_int32 :: error_flag
1614  sll_int32 :: error_flag2
1615  sll_real64, dimension(3, 3) :: p
1616  sll_real64, dimension(3, 1) :: bp1
1617  sll_real64, dimension(3, 1) :: bp2
1618  sll_int32, dimension(3) :: pivot
1619 
1620  ! Getting the indices of the vertices, have to compute first the center coos
1621  ccx = mesh%center_cartesian_coord(1, i_elmt)
1622  ccy = mesh%center_cartesian_coord(2, i_elmt)
1623  call sll_s_get_cell_vertices_index(ccx, ccy, mesh, e1, e2, e3)
1624 
1625  ! Getting the coordinates of the vertices of the element (in the hexmesh)
1626  x_v1 = mesh%cartesian_coord(1, e1)
1627  y_v1 = mesh%cartesian_coord(2, e1)
1628  x_v2 = mesh%cartesian_coord(1, e2)
1629  y_v2 = mesh%cartesian_coord(2, e2)
1630  x_v3 = mesh%cartesian_coord(1, e3)
1631  y_v3 = mesh%cartesian_coord(2, e3)
1632 
1633  ! Getting the coordinates of the mapped triangle
1634  call mesh%hex_to_aligned_pt(e1, transf, xp_v1, yp_v1)
1635  call mesh%hex_to_aligned_pt(e2, transf, xp_v2, yp_v2)
1636  call mesh%hex_to_aligned_pt(e3, transf, xp_v3, yp_v3)
1637 
1638  ! We now want to get A and B in AX + B = X', where X contains the values of
1639  ! the vertices in the hexmesh and X' the values of the mapped vertices,
1640  ! A = ( a11 a12, a21 a22) a 2x2 matrix et B=(b1 b2).
1641  ! For this we can solve two indepent linear systems of the form:
1642  ! P * A1 = Bp1, where the i-th line of 3x3 matrix P = (x_vi, y_vi, 1),
1643  ! A1 = (a11, a12, b1) and Bp1 = (xp_v1, xp_v2, xp_v3).
1644  ! The second linear system is P * A2 = Bp2. where P is the same matrix as
1645  ! in the previous system, A2 = (a21, a22, b2) and Bp2 is like Bp1 but
1646  ! using the y coordinates of the mapped points.
1647  ! We use LAPACK for the resolution.
1648 
1649  ! We first create the P matrices:
1650  p(1, 1:2) = (/x_v1, y_v1/)
1651  p(2, 1:2) = (/x_v2, y_v2/)
1652  p(3, 1:2) = (/x_v3, y_v3/)
1653  p(:, 3) = 1._f64
1654 
1655  ! Then the RHS vectors:
1656  bp1(1, 1) = xp_v1
1657  bp1(2, 1) = xp_v2
1658  bp1(3, 1) = xp_v3
1659  ! ... 2nd RHS vector:
1660  bp2(1, 1) = yp_v1
1661  bp2(2, 1) = yp_v2
1662  bp2(3, 1) = yp_v3
1663 
1664  ! Now we solve the system :
1665  error_flag = 0
1666  call dgesv(3, 1, p, 3, pivot, bp1, 3, error_flag)
1667  pivot(:) = 0
1668  p(1, 1:2) = (/x_v1, y_v1/)
1669  p(2, 1:2) = (/x_v2, y_v2/)
1670  p(3, 1:2) = (/x_v3, y_v3/)
1671  p(:, 3) = 1._f64
1672  call dgesv(3, 1, p, 3, pivot, bp2, 3, error_flag2)
1673  if ((error_flag .ne. 0) .or. (error_flag2 .ne. 0)) then
1674  print *, "LAPACK error flag for 1st sys = ", error_flag
1675  print *, "LAPACK error flag for 2nd sys = ", error_flag2
1676  sll_error('hex_to_aligned_elmt', "Error while calling DGESV from Lapack")
1677  end if
1678 
1679  transf_mata(1, 1) = bp1(1, 1)
1680  transf_mata(1, 2) = bp1(2, 1)
1681  transf_mata(2, 1) = bp2(1, 1)
1682  transf_mata(2, 2) = bp2(2, 1)
1683 
1684  transf_vecb(1) = bp1(3, 1)
1685  transf_vecb(2) = bp2(3, 1)
1686 
1687  end subroutine hex_to_aligned_elmt
1688 
1689  !---------------------------------------------------------------------------
1703  subroutine ref_to_hex_elmt(mesh, i_elmt, transf_matA, transf_vecB)
1704  class(sll_t_hex_mesh_2d), intent(in) :: mesh
1705  sll_int32, intent(in) :: i_elmt
1706  sll_real64, dimension(2, 2), intent(out) :: transf_mata
1707  sll_real64, dimension(2), intent(out) :: transf_vecb
1708  ! Local
1709  sll_int32 :: e1
1710  sll_int32 :: e2
1711  sll_int32 :: e3
1712  sll_int32 :: type
1713  sll_real64 :: x1, y1
1714  sll_real64 :: x_v1
1715  sll_real64 :: y_v1
1716 
1717  ! We get the cells center coordinates in order to get its vertices
1718  x1 = mesh%center_cartesian_coord(1, i_elmt)
1719  y1 = mesh%center_cartesian_coord(2, i_elmt)
1720  call sll_s_get_cell_vertices_index(x1, y1, mesh, e1, e2, e3)
1721 
1722  ! We get the lowest vertex coordinates (by default: e1's coordinates)
1723  x_v1 = mesh%cartesian_coord(1, e1)
1724  y_v1 = mesh%cartesian_coord(2, e1)
1725 
1726  ! Getting the cell type:
1727  call mesh%cell_type(i_elmt, type)
1728 
1729  ! We fill the matrices A and B:
1730  transf_vecb(1:2) = (/x_v1, y_v1/)
1731 
1732  if (type == 1) then
1733  transf_mata(1, 1) = 0.5_f64/real(mesh%num_cells, f64)
1734  transf_mata(1, 2) = -sll_p_sqrt3/2._f64/real(mesh%num_cells, f64)
1735  transf_mata(2, 1) = sll_p_sqrt3/2._f64/real(mesh%num_cells, f64)
1736  transf_mata(2, 2) = 0.5_f64/real(mesh%num_cells, f64)
1737  else
1738  transf_mata(1, 1) = 1._f64/real(mesh%num_cells, f64)
1739  transf_mata(1, 2) = 0._f64/real(mesh%num_cells, f64)
1740  transf_mata(2, 1) = 0._f64/real(mesh%num_cells, f64)
1741  transf_mata(2, 2) = 1._f64/real(mesh%num_cells, f64)
1742  end if
1743 
1744  end subroutine ref_to_hex_elmt
1745 
1746  !---------------------------------------------------------------------------
1757  subroutine ref_to_aligned_elmt(mesh, i_elmt, transf, transf_matA, transf_vecB)
1758  class(sll_t_hex_mesh_2d), intent(in) :: mesh
1759  sll_int32, intent(in) :: i_elmt
1760  character(len=*), intent(in) :: transf
1761  sll_real64, dimension(2, 2), intent(out) :: transf_mata
1762  sll_real64, dimension(2), intent(out) :: transf_vecb
1763  ! Local
1764  sll_real64, dimension(2, 2) :: transf_mata1
1765  sll_real64, dimension(2) :: transf_vecb1
1766  sll_real64, dimension(2, 2) :: transf_mata2
1767  sll_real64, dimension(2) :: transf_vecb2
1768 
1769  ! To compute the transformation ref->aligned we compute first the
1770  ! transformations ref->hex then hex->aligned:
1771  call mesh%ref_to_hex_elmt(i_elmt, transf_mata1, transf_vecb1)
1772  call mesh%hex_to_aligned_elmt(i_elmt, transf, transf_mata2, transf_vecb2)
1773 
1774  ! To get the composition of two linear combinations, we know that:
1775  ! A = A2*A1 B=A2*B1+B2
1776  transf_mata = matmul(transf_mata2, transf_mata1)
1777  transf_vecb = matmul(transf_mata2, transf_vecb1) + transf_vecb2
1778 
1779  end subroutine ref_to_aligned_elmt
1780 
1781  !-----------------------------------------------------------------------------
1785  subroutine sll_s_display_hex_mesh_2d(mesh)
1786  ! Displays mesh information on the terminal
1787  class(sll_t_hex_mesh_2d), intent(in) :: mesh
1788 
1789  write (*, "(/,(a))") '2D mesh : num_cells num_pts center_x1 &
1790  & center_x2 &
1791  & radius'
1792  write (*, "(10x,2(i6,9x),3(g13.3,1x))") mesh%num_cells, &
1793  mesh%num_pts_tot, &
1794  mesh%center_x1, &
1795  mesh%center_x2, &
1796  mesh%radius
1797  end subroutine sll_s_display_hex_mesh_2d
1798 
1805  subroutine sll_s_write_caid_files(mesh, transf, spline_deg)
1806  type(sll_t_hex_mesh_2d), pointer :: mesh
1807  character(len=*), intent(in) :: transf
1808  character(len=20), parameter :: name_nodes = "boxsplines_nodes.txt"
1809  character(len=23), parameter :: name_elemt = "boxsplines_elements.txt"
1810  character(len=24), parameter :: name_diri = "boxsplines_dirichlet.txt"
1811  sll_real64 :: x1, y1
1812  sll_real64 :: scale
1813  sll_int32 :: e1, e2, e3
1814  sll_int32 :: temp_e
1815  sll_int32 :: i, j
1816  sll_int32 :: k1, k2
1817  sll_int32 :: spline_deg
1818  sll_int32 :: nen
1819  sll_int32 :: num_pts_tot
1820  sll_int32 :: num_ele
1821  sll_int32 :: nei1, nei2, nei3
1822  sll_int32 :: num_cells_to_origin
1823  sll_int32 :: boundary
1824  sll_int32 :: dirichlet
1825  sll_int32 :: type
1826  sll_int32, parameter :: out_unit = 20
1827  sll_real64, dimension(2, 2) :: mata
1828  sll_real64, dimension(2) :: vecb
1829 
1830  ! Writing the nodes file....................
1831  open (unit=out_unit, file=name_nodes, action="write", status="replace")
1832 
1833  ! We first write the total number of points/nodes:
1834  num_pts_tot = mesh%num_pts_tot
1835  write (out_unit, "(i6)") num_pts_tot
1836  ! For every node...
1837  do i = 1, num_pts_tot
1838  k1 = mesh%hex_coord(1, i)
1839  k2 = mesh%hex_coord(2, i)
1840  num_cells_to_origin = sll_f_cells_to_origin(k1, k2)
1841  boundary = 0
1842  if (num_cells_to_origin .eq. mesh%num_cells) then
1843  boundary = 1
1844  end if
1845  ! Mapping to circle:
1846  if (transf .eq. "IDENTITY") then
1847  x1 = mesh%cartesian_coord(1, i)
1848  y1 = mesh%cartesian_coord(2, i)
1849  else
1850  call mesh%hex_to_aligned_pt(i, transf, x1, y1)
1851  end if
1852 
1853  !... we write the coordinates
1854  write (out_unit, "((i6),(a,1x),(g25.17),(a,1x),(g25.17))") boundary, &
1855  ",", &
1856  x1, &
1857  ",", &
1858  y1
1859  end do
1860 
1861  close (out_unit)
1862 
1863  ! Writing the elements file....................
1864  ! File containing general information about the cells and the transformation
1865  open (unit=out_unit, file=name_elemt, action="write", status="replace")
1866 
1867  ! We first write the total number of cells/elements:
1868  num_ele = mesh%num_triangles
1869  write (out_unit, "(i6)") num_ele
1870 
1871  ! The scale is fix here
1872  scale = 1._f64
1873  !... we write its global number
1874  write (out_unit, "(i6)") spline_deg
1875 
1876  ! For every element...
1877  do i = 1, num_ele
1878  !... we write its global number
1879  write (out_unit, "(i6)") sll_f_change_elements_notation(mesh, i)
1880  !... we write its type (1 or 2)
1881  call mesh%cell_type(i, type)
1882  write (out_unit, "(i6)") type
1883  !... we write the spline degree
1884  write (out_unit, "(i6)") spline_deg
1885  !... we write the scale of the element
1886  write (out_unit, "((f22.17),(a,1x))", advance='no') scale, ","
1887  !... we write its neighbours
1888  call mesh%get_neighbours(i, nei1, nei2, nei3)
1889  write (out_unit, "(3((i6),(a,1x)))", advance='no') &
1890  sll_f_change_elements_notation(mesh, nei1), ",", &
1891  sll_f_change_elements_notation(mesh, nei2), ",", &
1892  sll_f_change_elements_notation(mesh, nei3), ","
1893  !... we write the indices of the edges
1894  x1 = mesh%center_cartesian_coord(1, i)
1895  y1 = mesh%center_cartesian_coord(2, i)
1896  call sll_s_get_cell_vertices_index(x1, y1, mesh, e1, e2, e3)
1897  if (type .ne. 2) then
1898  temp_e = e2
1899  e2 = e3
1900  e3 = temp_e
1901  end if
1902  write (out_unit, "((i6),(a,1x),(i6),(a,1x),(i6))") e1, ",", e2, ",", e3
1903  !... we write the coordinate transformation (*)
1904  select case (transf)
1905  case ("IDENTITY")
1906  call mesh%ref_to_hex_elmt(i, mata, vecb)
1907  case ("TOKAMAK")
1908  call mesh%ref_to_aligned_elmt(i, transf, mata, vecb)
1909  case ("CIRCLE")
1910  call mesh%ref_to_aligned_elmt(i, transf, mata, vecb)
1911  case default
1912  print *, "ERROR in write_caid_files(): Not known transfromation."
1913  print *, " Options are : CIRCLE, TOKAMAK and IDENTITY."
1914  stop
1915  end select
1916 
1917  write (out_unit, "(5((f22.17), (a,1x)), (f22.17))") &
1918  mata(1, 1), ",", mata(1, 2), ",", mata(2, 1), ",", mata(2, 2), ",", &
1919  vecb(1), ",", vecb(2)
1920  end do
1921  print *, ""
1922  close (out_unit)
1923 
1924  ! (*) The coordinate transformation is the transformation from the reference
1925  ! element to the current cell. As the reference element is the 1st cell of
1926  ! an hexagonal mesh of radius 1, the transformation is only a rotation
1927  ! followed by a translation. Thus we only need 6 values to stock the
1928  ! transformation. 4 values for the matrix A and 2 for the vector v, where:
1929  ! Ax + b = x'. x being the reference coordinates and x' the coordinates of
1930  ! the current mesh.
1931  ! Reference coordinates: (0,0), (sqrt(3)/2, 0.5), (0,1)
1932 
1933  ! Writing the dirichlet file....................
1934  open (unit=out_unit, file=name_diri, action="write", status="replace")
1935 
1936  ! We first write the total number of cells/elements:
1937  num_ele = mesh%num_triangles
1938  write (out_unit, "(i6)") num_ele
1939 
1940  ! The number of splines non null in a cell depends on the degree
1941  nen = 3*spline_deg*spline_deg
1942  dirichlet = 1
1943 
1944  ! For every element...
1945  do i = 1, num_ele
1946  !... we write its global number
1947  write (out_unit, "(i6)") sll_f_change_elements_notation(mesh, i)
1948  !... we write the number of elements non-nul
1949  write (out_unit, "(i6)") nen
1950  do j = 1, nen
1951  !... we write 1 as at the moment they are all dirichlet
1952  write (out_unit, "((i6),(a,1x))", advance='no') dirichlet, ","
1953  end do
1954  write (out_unit, *) ""
1955  end do
1956  close (out_unit)
1957 
1958  end subroutine sll_s_write_caid_files
1959 
1960  !---------------------------------------------------------------------------
1967  subroutine write_hex_mesh_2d(mesh, name)
1968  class(sll_t_hex_mesh_2d), intent(in) :: mesh
1969  character(len=*), intent(in) :: name
1970 
1971  sll_int32 :: i
1972  sll_int32 :: num_pts_tot
1973  sll_int32 :: k1, k2
1974  sll_int32 :: out_unit
1975 
1976  open (file=name, status="replace", form="formatted", newunit=out_unit)
1977 
1978  num_pts_tot = mesh%num_pts_tot
1979 
1980  ! Optional writing every mesh point and its cartesian coordinates :
1981  ! write(*,"(/,(a))") 'hex mesh : num_pnt x1 x2'
1982 
1983  do i = 1, num_pts_tot
1984  k1 = mesh%global_to_hex1(i)
1985  k2 = mesh%global_to_hex2(i)
1986  write (out_unit, "(3(i6,1x),2(g13.3,1x))") i, &
1987  k1, &
1988  k2, &
1989  mesh%global_to_x1(i), &
1990  mesh%global_to_x2(i)
1991  end do
1992 
1993  close (out_unit)
1994  end subroutine write_hex_mesh_2d
1995 
1996  !---------------------------------------------------------------------------
2004  subroutine write_field_hex_mesh_xmf(mesh, field, name)
2005 
2006  class(sll_t_hex_mesh_2d), intent(in) :: mesh
2007  sll_real64, dimension(:), intent(in) :: field
2008  character(len=*), intent(in) :: name
2009 
2010  sll_int32 :: i
2011  sll_int32 :: num_triangles
2012  sll_int32 :: num_pts_tot
2013  sll_real64, allocatable :: coor(:, :)
2014  sll_int32, allocatable :: ntri(:, :)
2015  sll_int32 :: error
2016  sll_real64 :: x1, x2
2017 
2018  num_pts_tot = mesh%num_pts_tot
2019  num_triangles = mesh%num_triangles
2020  sll_allocate(coor(2, num_pts_tot), error)
2021  sll_allocate(ntri(3, num_triangles), error)
2022 
2023  do i = 1, num_pts_tot
2024  coor(1, i) = mesh%global_to_x1(i)
2025  coor(2, i) = mesh%global_to_x2(i)
2026  end do
2027 
2028  do i = 1, num_triangles
2029  x1 = mesh%center_cartesian_coord(1, i)
2030  x2 = mesh%center_cartesian_coord(2, i)
2031  call sll_s_get_cell_vertices_index( &
2032  x1, x2, mesh, &
2033  ntri(1, i), ntri(2, i), ntri(3, i))
2034  end do
2035 
2036  call sll_s_write_tri_mesh_xmf( &
2037  name, coor, ntri, &
2038  num_pts_tot, num_triangles, field, 'values')
2039 
2040  end subroutine write_field_hex_mesh_xmf
2041 
2042  !---------------------------------------------------------------------------
2050  subroutine write_hex_mesh_mtv(mesh, mtv_file)
2051  class(sll_t_hex_mesh_2d), intent(in) :: mesh
2052  character(len=*), intent(in) :: mtv_file
2053 
2054  sll_real64 :: coor(2, mesh%num_pts_tot)
2055  sll_int32 :: ntri(3, mesh%num_triangles)
2056  sll_real64 :: x1
2057  sll_real64 :: y1
2058  sll_int32 :: is1
2059  sll_int32 :: is2
2060  sll_int32 :: is3
2061  sll_int32 :: out_unit
2062  sll_int32 :: i
2063 
2064  open (file=mtv_file, status="replace", form="formatted", newunit=out_unit)
2065 
2066  !--- Drawing mesh ---
2067  write (out_unit, "(a)") "$DATA=CURVE3D"
2068  write (out_unit, "(a)") "%equalscale=T"
2069  write (out_unit, "(a)") "%toplabel='Mesh' "
2070 
2071  do i = 1, mesh%num_triangles
2072  x1 = mesh%center_cartesian_coord(1, i)
2073  y1 = mesh%center_cartesian_coord(2, i)
2074 
2075  call sll_s_get_cell_vertices_index(x1, y1, mesh, is1, is2, is3)
2076 
2077  ntri(1, i) = is1
2078  ntri(2, i) = is2
2079  ntri(3, i) = is3
2080 
2081  coor(1, is1) = mesh%global_to_x1(is1)
2082  coor(2, is1) = mesh%global_to_x2(is1)
2083  coor(1, is2) = mesh%global_to_x1(is2)
2084  coor(2, is2) = mesh%global_to_x2(is2)
2085  coor(1, is3) = mesh%global_to_x1(is3)
2086  coor(2, is3) = mesh%global_to_x2(is3)
2087 
2088  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2089  write (out_unit, "(3f10.5)") coor(:, ntri(2, i)), 0.
2090  write (out_unit, "(3f10.5)") coor(:, ntri(3, i)), 0.
2091  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2092  write (out_unit, *)
2093  end do
2094 
2095  !--- Number of knots and triangles
2096 
2097  write (out_unit, "(a)") "$DATA=CURVE3D"
2098  write (out_unit, "(a)") "%equalscale=T"
2099  write (out_unit, "(a)") "%toplabel='Numeber of knots and triangles' "
2100 
2101  do i = 1, mesh%num_triangles
2102  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2103  write (out_unit, "(3f10.5)") coor(:, ntri(2, i)), 0.
2104  write (out_unit, "(3f10.5)") coor(:, ntri(3, i)), 0.
2105  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2106  write (out_unit, *)
2107  end do
2108 
2109  do i = 1, mesh%num_triangles
2110  x1 = (coor(1, ntri(1, i)) &
2111  + coor(1, ntri(2, i)) &
2112  + coor(1, ntri(3, i)))/3.0_f64
2113  y1 = (coor(2, ntri(1, i)) &
2114  + coor(2, ntri(2, i)) &
2115  + coor(2, ntri(3, i)))/3.0_f64
2116  write (out_unit, "(a)", advance="no") "@text x1="
2117  write (out_unit, "(f8.5)", advance="no") x1
2118  write (out_unit, "(a)", advance="no") " y1="
2119  write (out_unit, "(f8.5)", advance="no") y1
2120  write (out_unit, "(a)", advance="no") " z1=0. lc=4 ll='"
2121  write (out_unit, "(i4)", advance="no") i
2122  write (out_unit, "(a)") "'"
2123  end do
2124 
2125  do i = 1, mesh%num_pts_tot
2126  x1 = coor(1, i)
2127  y1 = coor(2, i)
2128  write (out_unit, "(a)", advance="no") "@text x1="
2129  write (out_unit, "(g15.3)", advance="no") x1
2130  write (out_unit, "(a)", advance="no") " y1="
2131  write (out_unit, "(g15.3)", advance="no") y1
2132  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
2133  write (out_unit, "(i4)", advance="no") i
2134  write (out_unit, "(a)") "'"
2135  end do
2136 
2137  !--- Knots' number
2138  write (out_unit, *) "$DATA=CURVE3D"
2139  write (out_unit, *) "%equalscale=T"
2140  write (out_unit, *) "%toplabel='Knots number' "
2141 
2142  do i = 1, mesh%num_triangles
2143  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2144  write (out_unit, "(3f10.5)") coor(:, ntri(2, i)), 0.
2145  write (out_unit, "(3f10.5)") coor(:, ntri(3, i)), 0.
2146  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2147  write (out_unit, *)
2148  end do
2149 
2150  do i = 1, mesh%num_pts_tot
2151  x1 = coor(1, i)
2152  y1 = coor(2, i)
2153  write (out_unit, "(a)", advance="no") "@text x1="
2154  write (out_unit, "(g15.3)", advance="no") x1
2155  write (out_unit, "(a)", advance="no") " y1="
2156  write (out_unit, "(g15.3)", advance="no") y1
2157  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
2158  write (out_unit, "(i4)", advance="no") i
2159  write (out_unit, "(a)") "'"
2160  end do
2161 
2162  !--- Triangles' number
2163  write (out_unit, *) "$DATA=CURVE3D"
2164  write (out_unit, *) "%equalscale=T"
2165  write (out_unit, *) "%toplabel='Triangles number' "
2166 
2167  do i = 1, mesh%num_triangles
2168  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2169  write (out_unit, "(3f10.5)") coor(:, ntri(2, i)), 0.
2170  write (out_unit, "(3f10.5)") coor(:, ntri(3, i)), 0.
2171  write (out_unit, "(3f10.5)") coor(:, ntri(1, i)), 0.
2172  write (out_unit, *)
2173  end do
2174 
2175  do i = 1, mesh%num_triangles
2176  x1 = (coor(1, ntri(1, i)) &
2177  + coor(1, ntri(2, i)) &
2178  + coor(1, ntri(3, i)))/3.0_f64
2179  y1 = (coor(2, ntri(1, i)) &
2180  + coor(2, ntri(2, i)) &
2181  + coor(2, ntri(3, i)))/3.0_f64
2182  write (out_unit, "(a)", advance="no") "@text x1="
2183  write (out_unit, "(g15.3)", advance="no") x1
2184  write (out_unit, "(a)", advance="no") " y1="
2185  write (out_unit, "(g15.3)", advance="no") y1
2186  write (out_unit, "(a)", advance="no") " z1=0. lc=4 ll='"
2187  write (out_unit, "(i4)", advance="no") i
2188  write (out_unit, "(a)") "'"
2189  end do
2190 
2191  write (out_unit, *) "$END"
2192  close (out_unit)
2193 
2194  end subroutine write_hex_mesh_mtv
2195 
2196  !---------------------------------------------------------------------------
2200  subroutine sll_s_delete_hex_mesh_2d(mesh)
2201  class(sll_t_hex_mesh_2d), intent(inout) :: mesh
2202  sll_int32 :: ierr
2203 
2204  sll_deallocate(mesh%cartesian_coord, ierr)
2205  sll_deallocate(mesh%hex_coord, ierr)
2206  sll_deallocate(mesh%global_indices, ierr)
2207  if (mesh%EXTRA_TABLES .eq. 1) then
2208  sll_deallocate(mesh%center_cartesian_coord, ierr)
2209  sll_deallocate(mesh%center_index, ierr)
2210  sll_deallocate(mesh%edge_center_cartesian_coord, ierr)
2211  sll_deallocate(mesh%edge_center_index, ierr)
2212  end if
2213 
2214  end subroutine sll_s_delete_hex_mesh_2d
2215 
2216 #undef TEST_PRESENCE_AND_ASSIGN_VAL
2217 end module sll_m_hexagonal_meshes
sll_m_hexagonal_meshes::sll_o_delete
Definition: sll_m_hexagonal_meshes.F90:140
sll_m_constants
Fortran module where set some physical and mathematical constants.
Definition: sll_m_constants.F90:23
sll_m_constants::sll_p_sqrt3
real(kind=f64), parameter, public sll_p_sqrt3
Definition: sll_m_constants.F90:85
sll_m_constants::sll_p_epsilon_0
real(kind=f64), parameter, public sll_p_epsilon_0
Definition: sll_m_constants.F90:66
sll_m_hexagonal_meshes
Definition: sll_m_hexagonal_meshes.F90:23
sll_m_hexagonal_meshes::hex_to_aligned_pt
recursive subroutine hex_to_aligned_pt(mesh, ind, transf, x_new, y_new)
Computes the coordinate transformation hex->aligned for a point.
Definition: sll_m_hexagonal_meshes.F90:1517
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_delete_hex_mesh_2d
subroutine, public sll_s_delete_hex_mesh_2d(mesh)
Deletes an hexagonal mesh.
Definition: sll_m_hexagonal_meshes.F90:2201
sll_m_hexagonal_meshes::write_hex_mesh_mtv
subroutine write_hex_mesh_mtv(mesh, mtv_file)
Same as write_hex_mesh but output is mtv file.
Definition: sll_m_hexagonal_meshes.F90:2051
sll_m_hexagonal_meshes::global_to_hex1
function global_to_hex1(mesh, index)
Gives the 1st hexagonal coordinate to point of global index "index".
Definition: sll_m_hexagonal_meshes.F90:935
sll_m_hexagonal_meshes::sll_f_local_to_global
function, public sll_f_local_to_global(mesh, ref_index, local)
Transforms local index in a given reference local indexation to a global index from the mesh.
Definition: sll_m_hexagonal_meshes.F90:1064
sll_m_hexagonal_meshes::index_hex_to_global
subroutine index_hex_to_global(mesh, k1, k2, index_tab)
Finds the index in the global array of the point which coordinates where passed as parameter.
Definition: sll_m_hexagonal_meshes.F90:678
sll_m_hexagonal_meshes::eta2_cell_hex
real(kind=f64) function eta2_cell_hex(mesh, cell_num)
Computes the 2nd cartesian coordinate of the center of the cell.
Definition: sll_m_hexagonal_meshes.F90:770
sll_m_hexagonal_meshes::global_to_local
function global_to_local(mesh, ref_index, global)
Transforms global to local index in respect to a reference index.
Definition: sll_m_hexagonal_meshes.F90:1030
sll_m_hexagonal_meshes::hex_to_global
function hex_to_global(mesh, k1, k2)
Transform hexagonal coordinates to global index.
Definition: sll_m_hexagonal_meshes.F90:907
sll_m_hexagonal_meshes::eta2_cell_hex_two_arg
real(kind=f64) function eta2_cell_hex_two_arg(mesh, i, j)
NOT IMPLEMENTED.
Definition: sll_m_hexagonal_meshes.F90:803
sll_m_hexagonal_meshes::local_hex_to_global
function local_hex_to_global(mesh, k1_ref, k2_ref, local)
Same as sll_f_local_to_global but taking hexagonal coordinates as reference indexing.
Definition: sll_m_hexagonal_meshes.F90:1097
sll_m_hexagonal_meshes::get_neighbours
subroutine get_neighbours(mesh, cell_index, nei_1, nei_2, nei_3)
returns the indices of the neighbouring cells/triangles
Definition: sll_m_hexagonal_meshes.F90:1236
sll_m_hexagonal_meshes::write_hex_mesh_2d
subroutine write_hex_mesh_2d(mesh, name)
Writes the hexagonal mesh into a given file.
Definition: sll_m_hexagonal_meshes.F90:1968
sll_m_hexagonal_meshes::init_edge_center_triangle
subroutine init_edge_center_triangle(mesh)
<BRIEF_DESCRIPTION>
Definition: sll_m_hexagonal_meshes.F90:573
sll_m_hexagonal_meshes::hex_to_aligned_elmt
subroutine hex_to_aligned_elmt(mesh, i_elmt, transf, transf_matA, transf_vecB)
Computes the coordinate transformation hex->aligned for an element.
Definition: sll_m_hexagonal_meshes.F90:1597
sll_m_hexagonal_meshes::sll_f_change_elements_notation
function, public sll_f_change_elements_notation(mesh, i_elmt_old)
Function that allows to change from the current element notation to one more intuitive.
Definition: sll_m_hexagonal_meshes.F90:1371
sll_m_hexagonal_meshes::sll_f_cells_to_origin
integer(kind=i32) function, public sll_f_cells_to_origin(k1, k2)
Computes the number of cell from point to origin.
Definition: sll_m_hexagonal_meshes.F90:822
sll_m_hexagonal_meshes::sll_f_cart_to_hex2
function, public sll_f_cart_to_hex2(mesh, x1, x2)
Transform cartesian coordinates to 2ndst hexagonal coordinates.
Definition: sll_m_hexagonal_meshes.F90:1010
sll_m_hexagonal_meshes::ref_to_hex_elmt
subroutine ref_to_hex_elmt(mesh, i_elmt, transf_matA, transf_vecB)
Computes the coordinate transformation ref->hex for an element.
Definition: sll_m_hexagonal_meshes.F90:1704
sll_m_hexagonal_meshes::eta2_node_hex
real(kind=f64) function eta2_node_hex(mesh, i, j)
Computes the second cartesian coordinate of a given point.
Definition: sll_m_hexagonal_meshes.F90:737
sll_m_hexagonal_meshes::sll_s_hex_mesh_2d_init
subroutine, public sll_s_hex_mesh_2d_init(mesh, num_cells, radius, center_x1, center_x2, r1_x1, r1_x2, r2_x1, r2_x2, r3_x1, r3_x2, EXTRA_TABLES)
initializes a previously allocated 2d hex-mesh
Definition: sll_m_hexagonal_meshes.F90:250
sll_m_hexagonal_meshes::cell_type
subroutine cell_type(mesh, num_ele, val)
Computes the type of triangle of a given cell.
Definition: sll_m_hexagonal_meshes.F90:845
sll_m_hexagonal_meshes::global_to_x1
function global_to_x1(mesh, index)
Gives the 1st cartesian coordinate to point of global index "index".
Definition: sll_m_hexagonal_meshes.F90:963
sll_m_hexagonal_meshes::sll_s_write_caid_files
subroutine, public sll_s_write_caid_files(mesh, transf, spline_deg)
Writes files for CAID.
Definition: sll_m_hexagonal_meshes.F90:1806
sll_m_hexagonal_meshes::sll_f_cart_to_hex1
function, public sll_f_cart_to_hex1(mesh, x1, x2)
Transform cartesian coordinates to 1st hexagonal coordinates.
Definition: sll_m_hexagonal_meshes.F90:992
sll_m_hexagonal_meshes::eta1_cell_hex
function eta1_cell_hex(mesh, cell_num)
Computes the first cartesian coordinate of the center of the cell.
Definition: sll_m_hexagonal_meshes.F90:753
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_f_new_hex_mesh_2d
type(sll_t_hex_mesh_2d) function, pointer, public sll_f_new_hex_mesh_2d(num_cells, center_x1, center_x2, r11, r12, r21, r22, r31, r32, radius, EXTRA_TABLES)
Creates and initializes a new hexagonal mesh.
Definition: sll_m_hexagonal_meshes.F90:191
sll_m_hexagonal_meshes::ref_to_aligned_elmt
subroutine ref_to_aligned_elmt(mesh, i_elmt, transf, transf_matA, transf_vecB)
Computes the coordinate transformation ref->aligned for an element.
Definition: sll_m_hexagonal_meshes.F90:1758
sll_m_hexagonal_meshes::global_to_hex2
function global_to_hex2(mesh, index)
Gives the 2nd hexagonal coordinate to point of global index "index".
Definition: sll_m_hexagonal_meshes.F90:949
sll_m_hexagonal_meshes::eta1_node_hex
real(kind=f64) function eta1_node_hex(mesh, i, j)
Computes the first cartesian coordinate of a given point.
Definition: sll_m_hexagonal_meshes.F90:719
sll_m_hexagonal_meshes::init_center_points_triangle
subroutine init_center_points_triangle(mesh)
<BRIEF_DESCRIPTION>
Definition: sll_m_hexagonal_meshes.F90:479
sll_m_hexagonal_meshes::write_field_hex_mesh_xmf
subroutine write_field_hex_mesh_xmf(mesh, field, name)
Same as write_field_hex_mesh but output in xmf.
Definition: sll_m_hexagonal_meshes.F90:2005
sll_m_hexagonal_meshes::eta1_cell_hex_two_arg
real(kind=f64) function eta1_cell_hex_two_arg(mesh, i, j)
NOT IMPLEMENTED.
Definition: sll_m_hexagonal_meshes.F90:787
sll_m_hexagonal_meshes::global_to_x2
function global_to_x2(mesh, index)
Gives the 2nd cartesian coordinate to point of global index "index".
Definition: sll_m_hexagonal_meshes.F90:977
sll_m_hexagonal_meshes::sll_s_display_hex_mesh_2d
subroutine, public sll_s_display_hex_mesh_2d(mesh)
Displays hexagonal mesh in terminal.
Definition: sll_m_hexagonal_meshes.F90:1786
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_meshes_base
Definition: sll_m_meshes_base.F90:1
sll_m_tri_mesh_xmf
Definition: sll_m_tri_mesh_xmf.F90:1
sll_m_tri_mesh_xmf::sll_s_write_tri_mesh_xmf
subroutine, public sll_s_write_tri_mesh_xmf(filename, coor, ntri, nbs, nbt, field, label)
Definition: sll_m_tri_mesh_xmf.F90:18
sll_m_hexagonal_meshes::hex_mesh_2d_ptr
Definition: sll_m_hexagonal_meshes.F90:136
sll_m_hexagonal_meshes::sll_t_hex_mesh_2d
2d hexagonal mesh
Definition: sll_m_hexagonal_meshes.F90:63
sll_m_meshes_base::sll_c_mesh_2d_base
2D logical mesh
Definition: sll_m_meshes_base.F90:26
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1