Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_triangular_meshes.F90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------
2 ! SELALIB
3 !------------------------------------------------------------------------------
4 ! MODULE: sll_m_triangular_meshes
5 !
6 ! DESCRIPTION:
18 !------------------------------------------------------------------------------
19 module sll_m_triangular_meshes
20 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21 #include "sll_assert.h"
22 #include "sll_errors.h"
23 #include "sll_memory.h"
24 #include "sll_working_precision.h"
25 
26  use sll_m_boundary_condition_descriptors, only: &
27  sll_p_periodic
28 
29  use sll_m_constants, only: &
30  sll_p_pi
31 
32  use sll_m_hexagonal_meshes, only: &
33  sll_s_get_cell_vertices_index, &
34  sll_t_hex_mesh_2d
35 
36  implicit none
37 
38  public :: &
39  sll_s_analyze_triangular_mesh, &
40  sll_s_map_to_circle, &
41  sll_o_new_triangular_mesh_2d, &
42  sll_s_triangular_mesh_2d_init_from_file, &
43  sll_s_triangular_mesh_2d_init_from_square, &
44  sll_s_triangular_mesh_2d_free, &
45  sll_s_read_from_file, &
46  sll_o_create, &
47  sll_o_delete, &
48  sll_o_display, &
49  sll_t_triangular_mesh_2d, &
50  sll_s_write_triangular_mesh_mtv
51 
52  private
53 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
54 
56 ! vtaux - coordinate x vectors tan
57 ! vtauy - coordinate y vectors tan
58  type :: sll_t_triangular_mesh_2d
59 
60  sll_int32 :: num_nodes
61  sll_int32 :: num_triangles
62  sll_int32 :: num_edges
63  sll_int32 :: num_bound
64 
65  sll_real64, pointer :: coord(:, :)
66  sll_int32, pointer :: nodes(:, :)
67 
68  sll_real64 :: eta1_min
69  sll_real64 :: eta1_max
70  sll_real64 :: eta2_min
71  sll_real64 :: eta2_max
72 
73  sll_int32 :: nbcoti
74  sll_int32 :: nbtcot
75  sll_int32 :: nmxfr
76  sll_int32 :: nelfr
77  sll_int32 :: nmxsd
78  sll_int32 :: nctfrt
79  sll_real64 :: petitl
80  sll_real64 :: grandl
81 
82  sll_real64, dimension(:), pointer :: area
83  sll_int32, dimension(:), pointer :: refs
84  sll_int32, dimension(:), pointer :: reft
85  sll_int32, dimension(:, :), pointer :: nvois
86  sll_int32, dimension(:), pointer :: nusd
87  sll_int32, dimension(:), pointer :: npoel1
88  sll_int32, dimension(:), pointer :: npoel2
89  sll_int32, dimension(:), pointer :: krefro
90  sll_int32, dimension(:), pointer :: kctfro
91  sll_int32, dimension(:), pointer :: kelfro
92  sll_int32, dimension(:, :), pointer :: ksofro
93  sll_real64, dimension(:, :), pointer :: vnofro
94  sll_real64, dimension(:), pointer :: xmal1
95  sll_real64, dimension(:), pointer :: xmal2
96  sll_real64, dimension(:), pointer :: xmal3
97  sll_int32, dimension(:, :), pointer :: nuvac
98  sll_int32, dimension(:), pointer :: nugcv
99  sll_int32, dimension(:), pointer :: nbcov
100  sll_real64, dimension(:), pointer :: xlcod
101 
102  sll_real64, dimension(:), pointer :: vtaux
103  sll_real64, dimension(:), pointer :: vtauy
104  sll_int32, dimension(:), pointer :: nctfro
105  sll_int32, dimension(:), pointer :: nctfrp
106 
107  logical :: analyzed = .false.
108 
109  contains
110 
111  procedure, pass(mesh) :: global_to_x1
112  procedure, pass(mesh) :: global_to_x2
113 
114  end type sll_t_triangular_mesh_2d
115 
116  interface sll_o_create
117  module procedure initialize_triangular_mesh_2d
118  end interface sll_o_create
119 
120  interface sll_o_delete
121  module procedure sll_s_triangular_mesh_2d_free
122  end interface sll_o_delete
123 
124  interface sll_o_display
125  module procedure display_triangular_mesh_2d
126  end interface sll_o_display
127 
128  interface sll_o_new_triangular_mesh_2d
129  module procedure new_triangular_mesh_2d_from_file
130  module procedure new_triangular_mesh_2d_from_hex_mesh
131  module procedure new_triangular_mesh_2d_from_square
132  end interface sll_o_new_triangular_mesh_2d
133 
134 !Type: sll_t_triangular_mesh_2d
135 !
136 !nbs - nodes
137 !nbt - triangles
138 !nbtcot - edges
139 !nbcoti - inside edges
140 !num_bound - total references for BC
141 !nelfr - triangles on boundary
142 !coor - coords nodes
143 !area - areas elements
144 !refs - references nodes
145 !reft - references elements
146 !nodes - connectivity
147 !nvois - neighbors
148 !nvoif - number of edge in neighbor triangle
149 !nusd - references of subdomain
150 !petitl - smaller length
151 !grandl - bigger length
152 !ncfrt - edges on boundary
153 !nbcfli - edges inside with CFL problem
154 !nndfnt - nodes Dirichlet inside boubdary
155 !noefnt - nodes Dirichlet boundary
156 !irffnt - reference nodes Dirichlet
157 !nmxsd - sub domains number
158 
159 contains
160 
166  function new_triangular_mesh_2d_from_file(maafil) result(m)
167 
168  type(sll_t_triangular_mesh_2d), pointer :: m
169  character(len=*), intent(in) :: maafil
170 
171  sll_int32 :: ierr
172  sll_allocate(m, ierr)
173  call sll_s_read_from_file(m, maafil)
174 
175  call compute_areas(m)
176 
177  end function new_triangular_mesh_2d_from_file
178 
182  subroutine sll_s_triangular_mesh_2d_init_from_file(m, maafil)
183 
184  type(sll_t_triangular_mesh_2d) :: m
185  character(len=*), intent(in) :: maafil
186 
187  call sll_s_read_from_file(m, maafil)
188  call compute_areas(m)
189 
190  end subroutine sll_s_triangular_mesh_2d_init_from_file
191 
197  function new_triangular_mesh_2d_from_hex_mesh(hex_mesh) result(tri_mesh)
198 
199  type(sll_t_hex_mesh_2d), intent(in), pointer :: hex_mesh
200  type(sll_t_triangular_mesh_2d), pointer :: tri_mesh
201 
202  sll_int32 :: ierr
203  sll_int32 :: is1, iv1
204  sll_int32 :: is2, iv2
205  sll_int32 :: is3, iv3
206  sll_int32 :: i
207  sll_int32 :: nctfr
208  sll_real64 :: x1
209  sll_real64 :: y1
210  sll_real64 :: xa, xb, xc
211  sll_real64 :: ya, yb, yc
212  sll_real64 :: det
213 
214  sll_allocate(tri_mesh, ierr)
215 
216  tri_mesh%num_nodes = hex_mesh%num_pts_tot
217  tri_mesh%num_triangles = hex_mesh%num_triangles
218  tri_mesh%nmxfr = 1
219  tri_mesh%nmxsd = 1
220  sll_allocate(tri_mesh%coord(1:2, tri_mesh%num_nodes), ierr)
221  sll_allocate(tri_mesh%nodes(1:3, 1:tri_mesh%num_triangles), ierr)
222  sll_allocate(tri_mesh%refs(tri_mesh%num_nodes), ierr)
223  sll_allocate(tri_mesh%nvois(1:3, 1:tri_mesh%num_triangles), ierr)
224  tri_mesh%refs = 0
225  tri_mesh%nvois = 0
226 
227  do i = 1, hex_mesh%num_pts_tot
228  tri_mesh%coord(1, i) = hex_mesh%global_to_x1(i)
229  tri_mesh%coord(2, i) = hex_mesh%global_to_x2(i)
230  end do
231 
232  nctfr = 0
233  do i = 1, hex_mesh%num_triangles
234 
235  x1 = hex_mesh%center_cartesian_coord(1, i)
236  y1 = hex_mesh%center_cartesian_coord(2, i)
237 
238  call sll_s_get_cell_vertices_index(x1, y1, hex_mesh, is1, is2, is3)
239  call hex_mesh%get_neighbours(i, iv1, iv2, iv3)
240 
241  xa = tri_mesh%coord(1, is1)
242  ya = tri_mesh%coord(2, is1)
243  xb = tri_mesh%coord(1, is2)
244  yb = tri_mesh%coord(2, is2)
245  xc = tri_mesh%coord(1, is3)
246  yc = tri_mesh%coord(2, is3)
247 
248  det = 2.0_f64*((xb - xa)*(yc - ya) - (xc - xa)*(yb - ya))
249 
250  if (det > 0) then
251  tri_mesh%nodes(:, i) = [is1, is2, is3]
252  tri_mesh%nvois(:, i) = [iv1, iv2, iv3]
253  else
254  tri_mesh%nodes(:, i) = [is1, is3, is2]
255  tri_mesh%nvois(:, i) = [iv3, iv2, iv1]
256  end if
257 
258  if (tri_mesh%nvois(1, i) < 0) then
259  tri_mesh%refs(tri_mesh%nodes(1, i)) = 1
260  tri_mesh%refs(tri_mesh%nodes(2, i)) = 1
261  nctfr = nctfr + 1
262  end if
263  if (tri_mesh%nvois(2, i) < 0) then
264  tri_mesh%refs(tri_mesh%nodes(2, i)) = 1
265  tri_mesh%refs(tri_mesh%nodes(3, i)) = 1
266  nctfr = nctfr + 1
267  end if
268  if (tri_mesh%nvois(3, i) < 0) then
269  tri_mesh%refs(tri_mesh%nodes(3, i)) = 1
270  tri_mesh%refs(tri_mesh%nodes(1, i)) = 1
271  nctfr = nctfr + 1
272  end if
273 
274  end do
275 
276  call compute_areas(tri_mesh)
277 
278  end function new_triangular_mesh_2d_from_hex_mesh
279 
292  function new_triangular_mesh_2d_from_square(nc_eta1, &
293  eta1_min, &
294  eta1_max, &
295  nc_eta2, &
296  eta2_min, &
297  eta2_max, &
298  bc) result(m)
299 
300  type(sll_t_triangular_mesh_2d), pointer :: m
301  sll_int32, intent(in) :: nc_eta1
302  sll_real64, intent(in) :: eta1_min
303  sll_real64, intent(in) :: eta1_max
304  sll_int32, intent(in) :: nc_eta2
305  sll_real64, intent(in) :: eta2_min
306  sll_real64, intent(in) :: eta2_max
307  sll_int32, optional :: bc
308 
309  sll_int32 :: ierr
310  sll_allocate(m, ierr)
311 
312  if (present(bc)) then
313  sll_assert(bc == sll_p_periodic)
314  end if
315 
316  call initialize_triangular_mesh_2d(m, &
317  nc_eta1, &
318  eta1_min, &
319  eta1_max, &
320  nc_eta2, &
321  eta2_min, &
322  eta2_max)
323 
324  call compute_areas(m)
325 
326  end function new_triangular_mesh_2d_from_square
327 
337  subroutine sll_s_triangular_mesh_2d_init_from_square(m, &
338  nc_eta1, eta1_min, eta1_max, nc_eta2, &
339  eta2_min, eta2_max, bc)
340 
341  type(sll_t_triangular_mesh_2d) :: m
342  sll_int32, intent(in) :: nc_eta1
343  sll_real64, intent(in) :: eta1_min
344  sll_real64, intent(in) :: eta1_max
345  sll_int32, intent(in) :: nc_eta2
346  sll_real64, intent(in) :: eta2_min
347  sll_real64, intent(in) :: eta2_max
348  sll_int32, optional :: bc
349 
350  if (present(bc)) then
351  sll_assert(bc == sll_p_periodic)
352  end if
353 
354  call initialize_triangular_mesh_2d(m, &
355  nc_eta1, &
356  eta1_min, &
357  eta1_max, &
358  nc_eta2, &
359  eta2_min, &
360  eta2_max)
361 
362  call compute_areas(m)
363 
364  end subroutine sll_s_triangular_mesh_2d_init_from_square
365 
366  subroutine initialize_triangular_mesh_2d(mesh, &
367  nc_eta1, &
368  eta1_min, &
369  eta1_max, &
370  nc_eta2, &
371  eta2_min, &
372  eta2_max, &
373  bc)
374 
375  type(sll_t_triangular_mesh_2d) :: mesh
376  sll_int32, intent(in) :: nc_eta1
377  sll_real64, intent(in) :: eta1_min
378  sll_real64, intent(in) :: eta1_max
379  sll_int32, intent(in) :: nc_eta2
380  sll_real64, intent(in) :: eta2_min
381  sll_real64, intent(in) :: eta2_max
382  sll_int32, optional :: bc
383 
384 !*-----------------------------------------------------------------------
385 !* Mesh generation of rectangular domain
386 !* TRIANGLES P1 - Q4T.
387 !*-----------------------------------------------------------------------
388 
389  sll_int32 :: i, j
390  sll_int32 :: nd1, nd2, nd3
391  sll_int32 :: nelt
392  sll_int32 :: l, l1, ll1, ll2
393  sll_int32 :: in, iin, n, noeud, nsom, nxm, nym, neltot, ndd
394  sll_int32 :: nbox, nboy
395  sll_int32 :: iel, iev
396  sll_int32 :: nelin, nefro, nelfr, nhp, nlp
397  sll_int32 :: ierr
398 
399  sll_real64 :: xx1, xx2, pasx0, pasx1, pasy0, pasy1
400  sll_real64 :: alx, aly
401  sll_int32, allocatable :: nar(:)
402 
403  mesh%eta1_min = eta1_min
404  mesh%eta1_max = eta1_max
405  mesh%eta2_min = eta2_min
406  mesh%eta2_max = eta2_max
407 
408  nbox = nc_eta1 + 1
409  nboy = nc_eta2 + 1
410  ndd = max(nbox, nboy)
411  alx = eta1_max - eta1_min
412  aly = eta2_max - eta2_min
413 
414  nxm = nbox - 1 !number quadrangles x
415  nym = nboy - 1 !number quadrangles y
416  neltot = 2*nxm*nym !number total triangles
417  nsom = nbox*nboy !number nodes
418  noeud = nsom !number total nodes
419 
420  mesh%num_triangles = neltot
421  mesh%num_nodes = noeud
422  sll_clear_allocate(mesh%coord(1:2, 1:noeud), ierr)
423  sll_allocate(mesh%nodes(1:3, 1:neltot), ierr); mesh%nodes = 0
424  sll_allocate(mesh%nvois(3, 1:neltot), ierr); mesh%nvois = 0
425  sll_allocate(mesh%refs(1:noeud), ierr); mesh%refs = 0
426 
427  pasx0 = alx/real(nxm, f64)
428  pasx1 = pasx0*0.5_f64
429  pasy0 = aly/real(nym, f64)
430  pasy1 = pasy0*0.5_f64
431 
432  do n = 0, nym
433  in = nbox*n
434  do i = 1, nbox
435  iin = in + i
436  xx1 = eta1_min + alx - real(i - 1, f64)*pasx0
437  xx2 = eta2_min + real(n, f64)*pasy0
438  mesh%coord(1, iin) = xx1
439  mesh%coord(2, iin) = xx2
440  end do
441  end do
442 
443  nelt = 0
444  do l = 1, nym
445 
446  l1 = l - 1
447  ll1 = l1*nbox
448  ll2 = l*nbox
449 
450  do i = 1, nxm
451 
452  nd1 = ll1 + i
453  nd2 = ll2 + i
454  nd3 = ll2 + i + 1
455 
456  nelt = nelt + 1
457 
458  mesh%nodes(1, nelt) = nd1
459  mesh%nodes(2, nelt) = nd2
460  mesh%nodes(3, nelt) = nd1 + 1
461 
462  nelt = nelt + 1
463 
464  mesh%nodes(1, nelt) = nd1 + 1
465  mesh%nodes(2, nelt) = nd2
466  mesh%nodes(3, nelt) = nd3
467 
468  end do
469 
470  end do
471 
472  nefro = 4*(nbox - 2) + 4*(nboy - 2)
473  nelfr = 2*(nbox - 1) + 2*(nboy - 1) - 2
474  nelin = neltot - nelfr
475 
476  allocate (nar(2*(nbox + nboy)))
477 
478  nhp = nbox + 1
479  do i = 1, nbox
480  nar(i) = nhp - i
481  mesh%refs(nar(i)) = 1
482  end do
483 
484  nhp = nbox*nym
485  do i = 1, nbox
486  nar(i) = nhp + i
487  mesh%refs(nar(i)) = 3
488  end do
489 
490  nlp = nboy + 1
491  do i = 1, nboy
492  nar(i) = nbox*(nlp - i)
493  mesh%refs(nar(i)) = 4
494  end do
495 
496  do i = 1, nboy
497  nar(i) = i*nbox - nxm
498  mesh%refs(nar(i)) = 2
499  end do
500 
501  mesh%nmxfr = 0
502 
503  if (present(bc)) then
504 
505  do i = 1, nym
506  iel = (i - 1)*2*nym + 1
507  iev = iel + 2*nym - 1
508  mesh%nvois(1, iel) = iev
509  mesh%nvois(3, iev) = iel
510  end do
511 
512  do i = 1, nxm
513  iel = (i - 1)*2 + 1
514  iev = iel + 2*nxm*(nym - 1) + 1
515  mesh%nvois(3, iel) = iev
516  mesh%nvois(2, iev) = iel
517  end do
518 
519  nelfr = 0
520  nefro = 0
521 
522  else
523 
524  do i = 1, 2*(nbox - 1), 2
525  mesh%nvois(3, i) = -1
526  mesh%nvois(2, neltot - i + 1) = -3
527  end do
528 
529  mesh%nvois(1, 1) = -2
530  nhp = 2*(nbox - 1)
531  do i = 1, nboy - 2
532  j = i*nhp
533  mesh%nvois(3, j) = -4
534  mesh%nvois(1, j + 1) = -2
535  end do
536  mesh%nvois(3, neltot) = -4
537 
538  do i = 1, noeud
539  if (mesh%refs(i) > mesh%nmxfr) mesh%nmxfr = mesh%nmxfr + 1
540  end do
541 
542  end if
543 
544  deallocate (nar)
545 
546  end subroutine initialize_triangular_mesh_2d
547 
549  subroutine display_triangular_mesh_2d(mesh)
550  class(sll_t_triangular_mesh_2d), intent(in) :: mesh
551 
552  sll_real64 :: eta1_min, eta1_max, eta2_min, eta2_max
553  sll_int32 :: i, nctfrt, nbtcot
554 
555  eta1_min = minval(mesh%coord(1, :))
556  eta1_max = maxval(mesh%coord(1, :))
557  eta2_min = minval(mesh%coord(2, :))
558  eta2_max = maxval(mesh%coord(2, :))
559 
560  nctfrt = 0
561  do i = 1, mesh%num_triangles
562  if (mesh%nvois(1, i) < 0) nctfrt = nctfrt + 1
563  if (mesh%nvois(2, i) < 0) nctfrt = nctfrt + 1
564  if (mesh%nvois(3, i) < 0) nctfrt = nctfrt + 1
565  end do
566 
567  nbtcot = (3*mesh%num_triangles + nctfrt)/2
568 
569  write (*, "(/,(a))") '2D mesh : num_triangles num_nodes num_edges '
570  write (*, "(10x,3(i6,9x),/,'Frame',/,4(g13.3,1x))") &
571  mesh%num_triangles, &
572  mesh%num_nodes, &
573  nbtcot, &
574  eta1_min, &
575  eta1_max, &
576  eta2_min, &
577  eta2_max
578 
579  end subroutine display_triangular_mesh_2d
580 
581  function global_to_x1(mesh, i) result(x1)
582 ! Takes the global index of the point (see hex_to_global(...) for conventions)
583 ! returns the first coordinate (x1) on the cartesian basis
584  class(sll_t_triangular_mesh_2d) :: mesh
585  sll_int32 :: i
586  sll_real64 :: x1
587 
588  x1 = mesh%coord(1, i)
589 
590  end function global_to_x1
591 
592  function global_to_x2(mesh, i) result(x2)
593 ! Takes the global index of the point (see hex_to_global(...) for conventions)
594 ! returns the second coordinate (x2) on the cartesian basis
595  class(sll_t_triangular_mesh_2d) :: mesh
596  sll_int32 :: i
597  sll_real64 :: x2
598 
599  x2 = mesh%coord(2, i)
600 
601  end function global_to_x2
602 
603  subroutine sll_s_write_triangular_mesh_mtv(mesh, mtv_file)
604 
605  type(sll_t_triangular_mesh_2d), intent(in) :: mesh
606  character(len=*), intent(in) :: mtv_file
607 
608  sll_real64 :: x1, x2
609  sll_real64 :: y1, y2
610  sll_int32 :: out_unit
611  sll_int32 :: i, j, k, l
612 
613  open (file=mtv_file, status='replace', form='formatted', newunit=out_unit)
614 
615 !--- Trace du maillage ---
616 
617  write (out_unit, "(a)") "$DATA=CURVE3D"
618  write (out_unit, "(a)") "%equalscale=T"
619  write (out_unit, "(a)") "%toplabel='Maillage' "
620 
621  do i = 1, mesh%num_triangles
622 
623  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
624  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(2, i)), 0.
625  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(3, i)), 0.
626  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
627  write (out_unit, *)
628 
629  end do
630 
631 !--- Numeros des noeuds et des triangles
632 
633  write (out_unit, "(a)") "$DATA=CURVE3D"
634  write (out_unit, "(a)") "%equalscale=T"
635  write (out_unit, "(a)") "%toplabel='Numeros des noeuds et des triangles' "
636 
637  do i = 1, mesh%num_triangles
638  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
639  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(2, i)), 0.
640  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(3, i)), 0.
641  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
642  write (out_unit, *)
643  end do
644 
645  do i = 1, mesh%num_triangles
646  x1 = (mesh%coord(1, mesh%nodes(1, i)) &
647  + mesh%coord(1, mesh%nodes(2, i)) &
648  + mesh%coord(1, mesh%nodes(3, i)))/3.0_f64
649  y1 = (mesh%coord(2, mesh%nodes(1, i)) &
650  + mesh%coord(2, mesh%nodes(2, i)) &
651  + mesh%coord(2, mesh%nodes(3, i)))/3.0_f64
652  write (out_unit, "(a)", advance="no") "@text x1="
653  write (out_unit, "(f8.5)", advance="no") x1
654  write (out_unit, "(a)", advance="no") " y1="
655  write (out_unit, "(f8.5)", advance="no") y1
656  write (out_unit, "(a)", advance="no") " z1=0. lc=4 ll='"
657  write (out_unit, "(i4)", advance="no") i
658  write (out_unit, "(a)") "'"
659  end do
660 
661  do i = 1, mesh%num_nodes
662  x1 = mesh%coord(1, i)
663  y1 = mesh%coord(2, i)
664  write (out_unit, "(a)", advance="no") "@text x1="
665  write (out_unit, "(g15.3)", advance="no") x1
666  write (out_unit, "(a)", advance="no") " y1="
667  write (out_unit, "(g15.3)", advance="no") y1
668  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
669  write (out_unit, "(i4)", advance="no") i
670  write (out_unit, "(a)") "'"
671  end do
672 
673 !--- Numeros des noeuds
674 
675  write (out_unit, *) "$DATA=CURVE3D"
676  write (out_unit, *) "%equalscale=T"
677  write (out_unit, *) "%toplabel='Numeros des noeuds' "
678 
679  do i = 1, mesh%num_triangles
680  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
681  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(2, i)), 0.
682  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(3, i)), 0.
683  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
684  write (out_unit, *)
685  end do
686 
687  do i = 1, mesh%num_nodes
688  x1 = mesh%coord(1, i)
689  y1 = mesh%coord(2, i)
690  write (out_unit, "(a)", advance="no") "@text x1="
691  write (out_unit, "(g15.3)", advance="no") x1
692  write (out_unit, "(a)", advance="no") " y1="
693  write (out_unit, "(g15.3)", advance="no") y1
694  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
695  write (out_unit, "(i4)", advance="no") i
696  write (out_unit, "(a)") "'"
697  end do
698 
699 !--- Numeros des triangles
700 
701  write (out_unit, *) "$DATA=CURVE3D"
702  write (out_unit, *) "%equalscale=T"
703  write (out_unit, *) "%toplabel='Numeros des triangles' "
704 
705  do i = 1, mesh%num_triangles
706  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
707  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(2, i)), 0.
708  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(3, i)), 0.
709  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nodes(1, i)), 0.
710  write (out_unit, *)
711  end do
712 
713  do i = 1, mesh%num_triangles
714  x1 = (mesh%coord(1, mesh%nodes(1, i)) &
715  + mesh%coord(1, mesh%nodes(2, i)) &
716  + mesh%coord(1, mesh%nodes(3, i)))/3.0_f64
717  y1 = (mesh%coord(2, mesh%nodes(1, i)) &
718  + mesh%coord(2, mesh%nodes(2, i)) &
719  + mesh%coord(2, mesh%nodes(3, i)))/3.0_f64
720  write (out_unit, "(a)", advance="no") "@text x1="
721  write (out_unit, "(g15.3)", advance="no") x1
722  write (out_unit, "(a)", advance="no") " y1="
723  write (out_unit, "(g15.3)", advance="no") y1
724  write (out_unit, "(a)", advance="no") " z1=0. lc=4 ll='"
725  write (out_unit, "(i4)", advance="no") i
726  write (out_unit, "(a)") "'"
727  end do
728 
729  write (out_unit, *)
730 
731  write (out_unit, *) "$DATA=CURVE3D"
732  write (out_unit, *) "%equalscale=T"
733  write (out_unit, *) "%toplabel='References des frontieres' "
734 
735  do i = 1, mesh%num_triangles
736  write (out_unit, *) mesh%coord(1:2, mesh%nodes(1, i)), 0.0
737  write (out_unit, *) mesh%coord(1:2, mesh%nodes(2, i)), 0.0
738  write (out_unit, *) mesh%coord(1:2, mesh%nodes(3, i)), 0.0
739  write (out_unit, *) mesh%coord(1:2, mesh%nodes(1, i)), 0.0
740  write (out_unit, *)
741  end do
742 
743  do i = 1, mesh%num_triangles
744  do j = 1, 3
745  if (mesh%nvois(j, i) < 0) then
746  k = mod(j - 1, 3) + 1 !Premier sommet
747  l = mod(j, 3) + 1 !Deuxieme sommet
748  x1 = mesh%coord(1, mesh%nodes(k, i))
749  y1 = mesh%coord(2, mesh%nodes(k, i))
750  x2 = mesh%coord(1, mesh%nodes(l, i))
751  y2 = mesh%coord(2, mesh%nodes(l, i))
752  write (out_unit, "(a)", advance="no") "@text x1="
753  write (out_unit, "(g15.3)", advance="no") 0.5_f64*(x1 + x2)
754  write (out_unit, "(a)", advance="no") " y1="
755  write (out_unit, "(g15.3)", advance="no") 0.5_f64*(y1 + y2)
756  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
757  write (out_unit, "(i1)", advance="no") - mesh%nvois(j, i)
758  write (out_unit, "(a)") "'"
759  end if
760  end do
761  end do
762 
763  write (out_unit, *)
764 
765  write (out_unit, *) "$DATA=CURVE3D"
766  write (out_unit, *) "%equalscale=T"
767  write (out_unit, *) "%toplabel='References des noeuds' "
768 
769  do i = 1, mesh%num_triangles
770  write (out_unit, *) mesh%coord(1:2, mesh%nodes(1, i)), 0.0
771  write (out_unit, *) mesh%coord(1:2, mesh%nodes(2, i)), 0.0
772  write (out_unit, *) mesh%coord(1:2, mesh%nodes(3, i)), 0.0
773  write (out_unit, *) mesh%coord(1:2, mesh%nodes(1, i)), 0.0
774  write (out_unit, *)
775  end do
776 
777  do i = 1, mesh%num_nodes
778  x1 = mesh%coord(1, i)
779  y1 = mesh%coord(2, i)
780  write (out_unit, "(a)", advance="no") "@text x1="
781  write (out_unit, "(g15.3)", advance="no") x1
782  write (out_unit, "(a)", advance="no") " y1="
783  write (out_unit, "(g15.3)", advance="no") y1
784  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
785  write (out_unit, "(i1)", advance="no") mesh%refs(i)
786  write (out_unit, "(a)") "'"
787  end do
788 
789  if (mesh%analyzed) then
790 
791  write (out_unit, *)
792 
793  write (out_unit, *) "$DATA=CURVE3D"
794  write (out_unit, *) "%equalscale=T"
795  write (out_unit, *) "%toplabel='Polygones de Voronoi'"
796 
797  do i = 1, mesh%num_triangles
798  write (out_unit, *) "%linetype = 1 # Solid Linetype (default=1)"
799  write (out_unit, *) "%linewidth = 1 # Linewidth (default=1)"
800  write (out_unit, *) "%linecolor = 1 # Line Color (default=1)"
801  do j = 1, 3
802  if (mesh%nvois(j, i) > 0) then
803  call get_cell_center(mesh, i, x1, y1)
804  write (out_unit, *) x1, y1, 0.
805  call get_cell_center(mesh, mesh%nvois(j, i), x1, y1)
806  write (out_unit, *) x1, y1, 0.
807  write (out_unit, *)
808  end if
809  end do
810  end do
811 
812  do i = 1, mesh%num_triangles
813 
814  write (out_unit, *) "%linetype = 1 # Solid Linetype (default=1)"
815  write (out_unit, *) "%linewidth = 1 # Linewidth (default=1)"
816  write (out_unit, *) "%linecolor = 2 # Line Color (default=1)"
817  write (out_unit, *) mesh%coord(1:2, mesh%nodes(1, i)), 0.
818  write (out_unit, *) mesh%coord(1:2, mesh%nodes(2, i)), 0.
819  write (out_unit, *) mesh%coord(1:2, mesh%nodes(3, i)), 0.
820  write (out_unit, *) mesh%coord(1:2, mesh%nodes(1, i)), 0.
821  write (out_unit, *)
822 
823  end do
824 
825  write (out_unit, *)
826  write (out_unit, *) "$DATA=CURVE3D"
827  write (out_unit, *) "%equalscale=T"
828  write (out_unit, *) "%toplabel='Numeros des noeuds et des cotes' "
829 
830  do i = 1, mesh%nbtcot
831  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nuvac(1, i)), 0.
832  write (out_unit, "(3f10.5)") mesh%coord(:, mesh%nuvac(2, i)), 0.
833  write (out_unit, *)
834  x1 = 0.5_f64*(mesh%coord(1, mesh%nuvac(1, i)) + mesh%coord(1, mesh%nuvac(2, i)))
835  y1 = 0.5_f64*(mesh%coord(2, mesh%nuvac(1, i)) + mesh%coord(2, mesh%nuvac(2, i)))
836  write (out_unit, "(a)", advance="no") "@text x1="
837  write (out_unit, "(g15.3)", advance="no") x1
838  write (out_unit, "(a)", advance="no") " y1="
839  write (out_unit, "(g15.3)", advance="no") y1
840  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
841  write (out_unit, "(i4)", advance="no") i
842  write (out_unit, "(a)") "'"
843  end do
844 
845  do i = 1, mesh%num_nodes
846  x1 = mesh%coord(1, i)
847  y1 = mesh%coord(2, i)
848  write (out_unit, "(a)", advance="no") "@text x1="
849  write (out_unit, "(g15.3)", advance="no") x1
850  write (out_unit, "(a)", advance="no") " y1="
851  write (out_unit, "(g15.3)", advance="no") y1
852  write (out_unit, "(a)", advance="no") " z1=0. lc=5 ll='"
853  write (out_unit, "(i4)", advance="no") i
854  write (out_unit, "(a)") "'"
855  end do
856 
857  end if
858 
859  write (out_unit, *) "$END"
860  close (out_unit)
861 
862  end subroutine sll_s_write_triangular_mesh_mtv
863 
864  subroutine sll_s_triangular_mesh_2d_free(mesh)
865 
866  class(sll_t_triangular_mesh_2d), intent(inout) :: mesh
867 
868  sll_assert(mesh%num_nodes > 0)
869 
870  end subroutine sll_s_triangular_mesh_2d_free
871 
872  subroutine sll_s_read_from_file(mesh, maafil)
873 
874  type(sll_t_triangular_mesh_2d) :: mesh
875  character(len=*), intent(in) :: maafil
876 
877  sll_int32 :: nfmaa = 12
878  sll_int32 :: sll_err
879  sll_int32 :: i
880  sll_int32 :: j
881  integer :: nelin
882  integer :: nefro
883  integer :: nmaill
884  integer :: iout = 6
885  integer :: imxref
886 
887  write (iout, "(/////10x,'>>> Read mesh from file <<<'/)")
888  print *, " will open : ", maafil
889  open (nfmaa, file=maafil, status='OLD', err=80)
890  print *, " opened : ", maafil
891  write (*, 1050, advance='no') trim(maafil)
892 
893  write (iout, "(10x,'Open the file' &
894  & /10x,'Unit number',i3,' fichier ',a14)") &
895  & nfmaa, maafil
896 
897  read (nfmaa, *)
898  read (nfmaa, *) nmaill, imxref
899  read (nfmaa, *) mesh%num_nodes, &
900  mesh%num_triangles, &
901  mesh%nmxfr, &
902  mesh%nmxsd, &
903  nefro, &
904  nelin, &
905  mesh%nelfr
906 
907  sll_allocate(mesh%coord(1:2, mesh%num_nodes), sll_err)
908  sll_allocate(mesh%nodes(1:3, 1:mesh%num_triangles), sll_err)
909  sll_allocate(mesh%refs(mesh%num_nodes), sll_err)
910  sll_allocate(mesh%reft(mesh%num_triangles), sll_err)
911  sll_allocate(mesh%nusd(mesh%num_triangles), sll_err)
912  sll_allocate(mesh%nvois(3, mesh%num_triangles), sll_err)
913 
914  read (nfmaa, *) ((mesh%coord(i, j), i=1, 2), j=1, mesh%num_nodes)
915  read (nfmaa, *) (mesh%refs(i), i=1, mesh%num_nodes)
916  read (nfmaa, *) ((mesh%nodes(i, j), i=1, 3), j=1, mesh%num_triangles)
917  read (nfmaa, *) ((mesh%nvois(i, j), i=1, 3), j=1, mesh%num_triangles)
918  read (nfmaa, *) (mesh%nusd(i), i=1, mesh%num_triangles)
919 
920  close (nfmaa)
921 
922  write (iout, "(//,10x,'Nb of nodes : ',i10/ &
923  & ,10x,'Nb of triangles : ',i10/ &
924  & ,10x,'Nb max of BC : ',i10/ &
925  & ,10x,'Nb max of subdomain : ',i10/ &
926  & ,10x,/'Nb triangles with 1 node' &
927  & ,' on boundary : ',i10/)") mesh%num_nodes, &
928  mesh%num_triangles, mesh%nmxfr, mesh%nmxsd, nefro
929 
930  write (iout, "(//10x,'Nb elements inside : ',i10/ &
931  & ,10x,'Nb elements boundary : ',i10/)") nelin, mesh%nelfr
932 
933 1050 format(/' Read mesh from file ', a, ' ? Y')
934 
935  return
936 80 continue
937  sll_error('Read from file', 'Input file '//maafil//' not found.')
938 
939  end subroutine sll_s_read_from_file
940 
941  subroutine get_cell_center(mesh, iel, x1, x2)
942 
943  class(sll_t_triangular_mesh_2d) :: mesh
944  sll_int32, intent(in) :: iel
945  sll_real64, intent(out) :: x1
946  sll_real64, intent(out) :: x2
947 
948  sll_real64 :: xa, ya, xb, yb, xc, yc
949  sll_real64 :: det, syca, syba
950 
951  xa = mesh%coord(1, mesh%nodes(1, iel))
952  ya = mesh%coord(2, mesh%nodes(1, iel))
953  xb = mesh%coord(1, mesh%nodes(2, iel))
954  yb = mesh%coord(2, mesh%nodes(2, iel))
955  xc = mesh%coord(1, mesh%nodes(3, iel))
956  yc = mesh%coord(2, mesh%nodes(3, iel))
957 
958  det = 2.0_f64*((xb - xa)*(yc - ya) - (xc - xa)*(yb - ya))
959  syca = (yc - ya)*(xb*xb - xa*xa + yb*yb - ya*ya)
960  syba = (xb - xa)*(xc*xc - xa*xa + yc*yc - ya*ya)
961 
962  x1 = (syca - (yb - ya)*(xc*xc - xa*xa + yc*yc - ya*ya))/det
963  x2 = (syba - (xc - xa)*(xb*xb - xa*xa + yb*yb - ya*ya))/det
964 
965  end subroutine get_cell_center
966 
972  subroutine sll_s_map_to_circle(mesh, num_cells, order)
973 
974  class(sll_t_triangular_mesh_2d), intent(inout) :: mesh
975  sll_int32, intent(in) :: num_cells
976  sll_int32, optional :: order
977 
978  sll_int32 :: i, j, cell, layer
979  sll_real64 :: r, alpha
980 
981  if (present(order)) then
982  layer = order
983  else
984  layer = num_cells
985  end if
986 
987  i = 2
988  do cell = 1, num_cells
989  if (cell > num_cells - layer) then
990  r = real(cell, f64)*1.0_f64/real(num_cells, f64)
991  alpha = sll_p_pi/6.0_f64
992  do j = 1, cell*6
993  mesh%coord(1, i + j - 1) = r*cos(alpha)
994  mesh%coord(2, i + j - 1) = r*sin(alpha)
995  alpha = alpha + sll_p_pi/real(3*cell, f64)
996  end do
997  end if
998  i = i + cell*6
999  end do
1000 
1001  end subroutine sll_s_map_to_circle
1002 
1003 !=======================================================================
1004 
1005  subroutine compute_areas(mesh)
1006 
1007  type(sll_t_triangular_mesh_2d), intent(inout) :: mesh
1008 
1009  integer, dimension(:), allocatable :: indc
1010 
1011  sll_int32 :: i, j, k
1012  sll_int32 :: it, is, is1, is2, is3
1013  sll_int32 :: ntmp, id1, nct, iel, ind, iel1, nel
1014  sll_int32 :: i1, i2, i3
1015  sll_int32 :: jel1, jel2, jel3, nel1, nel2, nel3
1016 
1017  real(8) :: airtot
1018  real(8) :: xlml, xlmu, ylml, ylmu
1019  real(8) :: lx1, lx2, ly1, ly2
1020 
1021  xlml = minval(mesh%coord(1, :))
1022  xlmu = maxval(mesh%coord(1, :))
1023  ylml = minval(mesh%coord(2, :))
1024  ylmu = maxval(mesh%coord(2, :))
1025 
1026  mesh%petitl = 1.d-04*min(xlmu - xlml, ylmu - ylml)/sqrt(real(mesh%num_nodes, 8))
1027  mesh%grandl = 1.d+04*max(xlmu - xlml, ylmu - ylmu)
1028 
1029  allocate (mesh%area(mesh%num_triangles)); mesh%area = 0.0_f64
1030 
1031  airtot = 0._f64
1032 
1033  do it = 1, mesh%num_triangles
1034 
1035  lx1 = mesh%coord(1, mesh%nodes(2, it)) - mesh%coord(1, mesh%nodes(1, it))
1036  ly1 = mesh%coord(2, mesh%nodes(3, it)) - mesh%coord(2, mesh%nodes(1, it))
1037  lx2 = mesh%coord(1, mesh%nodes(3, it)) - mesh%coord(1, mesh%nodes(1, it))
1038  ly2 = mesh%coord(2, mesh%nodes(2, it)) - mesh%coord(2, mesh%nodes(1, it))
1039 
1040  mesh%area(it) = 0.5_f64*abs(lx1*ly1 - lx2*ly2)
1041 
1042  if (mesh%area(it) <= 0._f64) then
1043  write (6, *) " Triangle : ", it
1044  write (6, *) mesh%nodes(1, it), ":", mesh%coord(1:2, mesh%nodes(1, it))
1045  write (6, *) mesh%nodes(2, it), ":", mesh%coord(1:2, mesh%nodes(2, it))
1046  write (6, *) mesh%nodes(3, it), ":", mesh%coord(1:2, mesh%nodes(3, it))
1047  stop "triangle negative area"
1048  end if
1049 
1050  airtot = airtot + mesh%area(it)
1051 
1052  end do
1053 
1054 ! --- triangles with one node in common -----------
1055 
1056 !search for elements with a common summit
1057 !creation of table npoel1(i+1) containing the number of
1058 !triangles having node i in common
1059 
1060  allocate (mesh%npoel1(mesh%num_nodes + 1))
1061 
1062  mesh%npoel1 = 0
1063  do i = 1, mesh%num_triangles
1064  is1 = mesh%nodes(1, i)
1065  is2 = mesh%nodes(2, i)
1066  is3 = mesh%nodes(3, i)
1067  mesh%npoel1(is1 + 1) = mesh%npoel1(is1 + 1) + 1
1068  mesh%npoel1(is2 + 1) = mesh%npoel1(is2 + 1) + 1
1069  mesh%npoel1(is3 + 1) = mesh%npoel1(is3 + 1) + 1
1070  end do
1071 
1072 ! table npoel1 becomes the table giving the address
1073 ! in npoel2 of the last element in the sequence of triangles
1074 ! common to a node
1075 
1076  mesh%npoel1(1) = 0
1077  do i = 3, mesh%num_nodes + 1
1078  mesh%npoel1(i) = mesh%npoel1(i - 1) + mesh%npoel1(i)
1079  end do
1080 
1081 ! creation of the npoel2 table containing sequentially the
1082 ! triangles numbers with a common node
1083 ! the first triangle leaning on node i is
1084 ! address by "npoel1(i)+1"
1085 ! the number of triangles having node i in common is
1086 ! "npoel1(i+1)-npoel1(i)"
1087 
1088  allocate (mesh%npoel2(mesh%npoel1(mesh%num_nodes + 1)))
1089  allocate (indc(mesh%num_nodes))
1090 
1091  indc = 1
1092 
1093  do it = 1, mesh%num_triangles
1094  do k = 1, 3
1095  is = mesh%nodes(k, it)
1096  mesh%npoel2(mesh%npoel1(is) + indc(is)) = it
1097  indc(is) = indc(is) + 1
1098  end do
1099  end do
1100 
1101  deallocate (indc)
1102 
1103 ! --- search neighbors
1104 
1105  do iel = 1, mesh%num_triangles
1106 
1107  is1 = mesh%nodes(1, iel)
1108  is2 = mesh%nodes(2, iel)
1109  is3 = mesh%nodes(3, iel)
1110 
1111  !nested loops on elements pointing towards
1112  !the 2 end nodes of the considered ridge
1113  !The neighbour is the common triangle (except iel)
1114  !first ridge (between vertex is1 and is2)
1115 
1116  nel1 = mesh%npoel1(is1 + 1) - mesh%npoel1(is1) !nb de triangles communs a is1
1117  nel2 = mesh%npoel1(is2 + 1) - mesh%npoel1(is2) !nb de triangles communs a is2
1118 
1119  loop1: do i1 = 1, nel1
1120  jel1 = mesh%npoel2(mesh%npoel1(is1) + i1) !premier triangle is1
1121  if (jel1 .ne. iel) then
1122  do i2 = 1, nel2
1123  jel2 = mesh%npoel2(mesh%npoel1(is2) + i2)
1124  if (jel2 == jel1) then
1125  mesh%nvois(1, iel) = jel1
1126  exit loop1
1127  end if
1128  end do
1129  end if
1130  end do loop1
1131 
1132  ! ... second edge (is2-is3)
1133 
1134  nel2 = mesh%npoel1(is2 + 1) - mesh%npoel1(is2)
1135  nel3 = mesh%npoel1(is3 + 1) - mesh%npoel1(is3)
1136 
1137  loop2: do i2 = 1, nel2
1138  jel2 = mesh%npoel2(mesh%npoel1(is2) + i2)
1139  if (jel2 /= iel) then
1140  do i3 = 1, nel3
1141  jel3 = mesh%npoel2(mesh%npoel1(is3) + i3)
1142  if (jel3 == jel2) then
1143  mesh%nvois(2, iel) = jel2
1144  exit loop2
1145  end if
1146  end do
1147  end if
1148  end do loop2
1149 
1150  ! ... third edge (is3-is1)
1151 
1152  nel3 = mesh%npoel1(is3 + 1) - mesh%npoel1(is3)
1153  nel1 = mesh%npoel1(is1 + 1) - mesh%npoel1(is1)
1154 
1155  loop3: do i3 = 1, nel3
1156  jel3 = mesh%npoel2(mesh%npoel1(is3) + i3)
1157  if (jel3 /= iel) then
1158  do i1 = 1, nel1
1159  jel1 = mesh%npoel2(mesh%npoel1(is1) + i1)
1160  if (jel1 == jel3) then
1161  mesh%nvois(3, iel) = jel3
1162  exit loop3
1163  end if
1164  end do
1165  end if
1166  end do loop3
1167 
1168  end do
1169 
1170 ! --- npoel2 ---------
1171 
1172  do is = 1, mesh%num_nodes
1173 
1174  nel = mesh%npoel1(is + 1) - mesh%npoel1(is)
1175 
1176  if (nel > 1) then
1177 
1178  !*** internal nodes (ref=0) ***
1179 
1180  if (mesh%refs(is) == 0) then
1181 
1182  ind = 1
1183  iel1 = mesh%npoel2(mesh%npoel1(is) + 1)
1184 
1185  loop4: do iel = 2, nel - 1
1186  do j = 1, 3
1187  if (mesh%nodes(j, iel1) == is) nct = mod(j + 1, 3) + 1
1188  end do
1189 
1190  iel1 = mesh%nvois(nct, iel1)
1191  do id1 = ind + 1, nel
1192  if (iel1 == mesh%npoel2(mesh%npoel1(is) + id1)) then
1193  ind = ind + 1
1194  ntmp = mesh%npoel2(mesh%npoel1(is) + ind)
1195  mesh%npoel2(mesh%npoel1(is) + ind) = iel1
1196  mesh%npoel2(mesh%npoel1(is) + id1) = ntmp
1197  cycle loop4
1198  end if
1199  end do
1200  end do loop4
1201 
1202  ! boundary nodes
1203 
1204  else
1205 
1206  ! --> first triangle
1207  loop5: do id1 = 1, nel
1208  iel1 = mesh%npoel2(mesh%npoel1(is) + id1)
1209  do j = 1, 3
1210  if (mesh%nvois(j, iel1) .le. 0 .and. mesh%nodes(j, iel1) == is) then
1211  ntmp = mesh%npoel2(mesh%npoel1(is) + 1)
1212  mesh%npoel2(mesh%npoel1(is) + 1) = iel1
1213  mesh%npoel2(mesh%npoel1(is) + id1) = ntmp
1214  exit loop5
1215  end if
1216  end do
1217  end do loop5
1218 
1219  if (nel > 2) then
1220  ind = 1
1221  iel1 = mesh%npoel2(mesh%npoel1(is) + 1)
1222 
1223  loop6: do iel = 2, nel - 1
1224  do j = 1, 3
1225  if (mesh%nodes(j, iel1) == is) then
1226  nct = mod(j + 1, 3) + 1
1227  end if
1228  end do
1229 
1230  iel1 = mesh%nvois(nct, iel1)
1231 
1232  do id1 = ind + 1, nel
1233  if (iel1 == mesh%npoel2(mesh%npoel1(is) + id1)) then
1234  ind = ind + 1
1235  ntmp = mesh%npoel2(mesh%npoel1(is) + ind)
1236  mesh%npoel2(mesh%npoel1(is) + ind) = iel1
1237  mesh%npoel2(mesh%npoel1(is) + id1) = ntmp
1238  cycle loop6
1239  end if
1240  end do
1241 
1242  end do loop6
1243 
1244  end if
1245 
1246  end if
1247 
1248  end if
1249 
1250  end do
1251 
1252  end subroutine compute_areas
1253 
1254 !**************************************************************
1255 
1256 !Subroutine: poclis
1257 !
1258 ! Calculation of the smoothing matrices associated with each
1259 ! mesh node.
1260 ! Required to calculate the components of E1,E2
1261 ! from the tangeantial components of E
1262 ! known on the dimensions of the triangles.
1263 ! Calculation of the components of the tangling unit vectors.
1264 !
1265 ! Input variables:
1266 !
1267 ! nuvac - PV numbers associated with the odds
1268 ! coor - coordinates of Delaunay nodes
1269 ! xlcod - Delaunay dimension length
1270 ! npoel1 - pointer of the npoel2 array
1271 ! npoel2 - number of triangles surrounding a node
1272 ! xmal1 - sum of the rates*rates surrounding a node (/det)
1273 ! xmal2 - sum of tauy*tauy surrounding a node (/det)
1274 ! xmal3 - sum of rates*tauy surrounding a node (/det)
1275 
1276  subroutine poclis(mesh, ncotcu, nuctfr)
1277 
1278  type(sll_t_triangular_mesh_2d) :: mesh
1279  integer, dimension(:) :: ncotcu
1280  integer, dimension(:) :: nuctfr
1281  logical :: lerr
1282  double precision :: det, s1, s2, s3, x21, y21
1283  double precision :: xa, ya, xb, yb
1284  integer :: nel1, nel2, indv1, indv2
1285  integer :: ic, nuctf, ind, nm1, nm2
1286  integer :: indn1, indn2, num1, num2
1287  integer :: n1, n2, ivois, nc, iel, nucti
1288  integer :: nbti, is
1289  integer, parameter :: iout = 6
1290  integer :: iac, nbc
1291 
1292 !======================================================================
1293 ! --- 1.0 --- Pointer of edge nodes -----
1294 
1295  do is = 1, mesh%num_nodes
1296  nbti = mesh%npoel1(is + 1) - mesh%npoel1(is)
1297  mesh%nbcov(is + 1) = nbti
1298  if (mesh%refs(is) /= 0) then
1299  mesh%nbcov(is + 1) = nbti + 1
1300  end if
1301  end do
1302 
1303  mesh%nbcov(1) = 0
1304  do is = 1, mesh%num_nodes
1305  mesh%nbcov(is + 1) = mesh%nbcov(is) + mesh%nbcov(is + 1)
1306  end do
1307 
1308  do is = 1, mesh%num_nodes
1309  nbc = mesh%nbcov(is + 1) - mesh%nbcov(is)
1310  iac = mesh%nbcov(is)
1311  end do
1312 
1313  nucti = 0
1314  nuctfr = 0
1315 
1316 ! ======================================================================
1317 ! --- 2.0 --- set edge numbers -----------------------------------
1318 
1319  do iel = 1, mesh%num_triangles
1320 
1321  do nc = 1, 3
1322 
1323  ivois = mesh%nvois(nc, iel)
1324  n1 = nc
1325  n2 = mod(nc, 3) + 1
1326 
1327  num1 = mesh%nodes(n1, iel)
1328  num2 = mesh%nodes(n2, iel)
1329  indn1 = mesh%npoel1(num1)
1330  indn2 = mesh%npoel1(num2)
1331  indv1 = mesh%nbcov(num1)
1332  indv2 = mesh%nbcov(num2)
1333  nel1 = mesh%npoel1(num1 + 1) - mesh%npoel1(num1)
1334  nel2 = mesh%npoel1(num2 + 1) - mesh%npoel1(num2)
1335 
1336  !internal edges
1337 
1338  if (ivois > iel) then
1339 
1340  nucti = nucti + 1
1341 
1342  !Numbers of edges with the same node
1343 
1344  do nm1 = 1, nel1
1345  if (mesh%npoel2(indn1 + nm1) == ivois) then
1346  mesh%nugcv(indv1 + nm1) = nucti
1347  end if
1348  end do
1349 
1350  do nm2 = 1, nel2
1351  if (mesh%npoel2(indn2 + nm2) == iel) then
1352  mesh%nugcv(indv2 + nm2) = nucti
1353  end if
1354  end do
1355 
1356  !number of triangles
1357 
1358  mesh%nuvac(1, nucti) = num1
1359  mesh%nuvac(2, nucti) = num2
1360 
1361  else if (ivois < 0) then !boundaries
1362 
1363  ind = -ivois
1364  nuctfr(ind) = nuctfr(ind) + 1
1365  nuctf = ncotcu(ind + 1) + nuctfr(ind)
1366  mesh%nugcv(indv1 + nel1 + 1) = nuctf
1367  mesh%nugcv(indv2 + nel2) = nuctf
1368  mesh%nuvac(1, nuctf) = num1
1369  mesh%nuvac(2, nuctf) = num2
1370 
1371  end if
1372 
1373  end do
1374 
1375  end do
1376 
1377 !======================================================================
1378 !----------- length triangles edges ------------------------
1379 
1380  do ic = 1, mesh%nbtcot
1381  xa = mesh%coord(1, mesh%nuvac(1, ic))
1382  ya = mesh%coord(2, mesh%nuvac(1, ic))
1383  xb = mesh%coord(1, mesh%nuvac(2, ic))
1384  yb = mesh%coord(2, mesh%nuvac(2, ic))
1385  mesh%xlcod(ic) = sqrt((xa - xb)*(xa - xb) + (ya - yb)*(ya - yb))
1386  end do
1387 
1388 !======================================================================
1389 !--- 4.0 --- matrices to compute electric field from potential --------
1390 
1391  do ic = 1, mesh%nbtcot
1392  n1 = mesh%nuvac(1, ic)
1393  n2 = mesh%nuvac(2, ic)
1394  x21 = (mesh%coord(1, n2) - mesh%coord(1, n1))/mesh%xlcod(ic)
1395  y21 = (mesh%coord(2, n2) - mesh%coord(2, n1))/mesh%xlcod(ic)
1396  s1 = x21*x21
1397  s2 = y21*y21
1398  s3 = x21*y21
1399  mesh%vtaux(ic) = x21
1400  mesh%vtauy(ic) = y21
1401  mesh%xmal1(n1) = mesh%xmal1(n1) + s1
1402  mesh%xmal2(n1) = mesh%xmal2(n1) + s2
1403  mesh%xmal3(n1) = mesh%xmal3(n1) + s3
1404  mesh%xmal1(n2) = mesh%xmal1(n2) + s1
1405  mesh%xmal2(n2) = mesh%xmal2(n2) + s2
1406  mesh%xmal3(n2) = mesh%xmal3(n2) + s3
1407  end do
1408 
1409 !--- 4.5 --- Normalization ---------------
1410  lerr = .false.
1411 
1412  do is = 1, mesh%num_nodes
1413  det = mesh%xmal1(is)*mesh%xmal2(is) - mesh%xmal3(is)*mesh%xmal3(is)
1414  if (det /= 0) then
1415  mesh%xmal1(is) = mesh%xmal1(is)/det
1416  mesh%xmal2(is) = mesh%xmal2(is)/det
1417  mesh%xmal3(is) = mesh%xmal3(is)/det
1418  else
1419  lerr = .true.
1420  end if
1421  end do
1422 
1423  if (lerr) then
1424  write (iout, 900)
1425  stop "poclis"
1426  end if
1427 
1428  if (lerr) then
1429  write (iout, 901)
1430  stop "poclis"
1431  end if
1432 
1433 900 format(//10x, 'Determinant des coefficients des matrices' &
1434  /10x, 'de lissage nul' &
1435  /10x, 'Il faut modifier le maillage'//)
1436 901 format(//10x, 'Le nombre de triangles communs a un noeud' &
1437  /10x, 'est superieur a 12' &
1438  /10x, 'Modifiez legerement le maillage'//)
1439 
1440  end subroutine poclis
1441 
1442 !========================================================================
1465  subroutine sll_s_analyze_triangular_mesh(mesh)
1466 
1467  integer, parameter :: iout = 6
1468  integer :: i, j, ifr
1469  type(sll_t_triangular_mesh_2d), intent(inout) :: mesh
1470 
1471  integer, dimension(:), allocatable :: nuctfr
1472  integer, dimension(:), allocatable :: ncotcu
1473 
1474  integer :: iel, is1, is2
1475  integer :: ict, ictcl, iref, jref
1476  integer :: keltmp, kcttmp, kretmp, ks1tmp, ks2tmp
1477  integer :: nbcot, ict1, ict2, ictcl1, ictcl2
1478  real(8) :: x1, y1, x2, y2
1479 
1480 !-----
1481 
1482  write (iout, "(//10x,'>>> Compute come quantities from mesh <<<'/)")
1483 
1484 ! --- Definition de nctfrt: le nombre de cotes frontieres
1485 
1486  mesh%nctfrt = 0
1487  do i = 1, mesh%num_triangles
1488  if (mesh%nvois(1, i) < 0) mesh%nctfrt = mesh%nctfrt + 1
1489  if (mesh%nvois(2, i) < 0) mesh%nctfrt = mesh%nctfrt + 1
1490  if (mesh%nvois(3, i) < 0) mesh%nctfrt = mesh%nctfrt + 1
1491  end do
1492 
1493 !*** total edges and internal (nbcoti,nbtcot)
1494 
1495  mesh%nbtcot = (3*mesh%num_triangles + mesh%nctfrt)/2
1496 
1497  mesh%nbcoti = mesh%nbtcot - mesh%nctfrt
1498 
1499 ! cumsum of edges ....
1500 
1501  allocate (ncotcu(mesh%nbcoti + mesh%nmxfr*mesh%nctfrt))
1502 
1503  ncotcu = 0
1504  ncotcu(1) = mesh%nbcoti
1505  ncotcu(2) = mesh%nbcoti
1506 
1507  do ifr = 1, mesh%nmxfr
1508  do i = 1, mesh%num_triangles
1509  do j = 1, 3
1510  if (mesh%nvois(j, i) == -ifr) then
1511  ncotcu(ifr + 2) = ncotcu(ifr + 2) + 1
1512  end if
1513  end do
1514  end do
1515  end do
1516 
1517  do ifr = 1, mesh%nmxfr
1518  ncotcu(ifr + 2) = ncotcu(ifr + 1) + ncotcu(ifr + 2)
1519  end do
1520 
1521  allocate (mesh%vtaux(mesh%nbtcot), mesh%vtauy(mesh%nbtcot))
1522  mesh%vtaux = 0._f64; mesh%vtauy = 0._f64
1523 
1524  allocate (mesh%nuvac(2, mesh%nbtcot)); mesh%nuvac = 0
1525 
1526  allocate (mesh%xlcod(mesh%nbtcot)); mesh%xlcod = 0.0_f64
1527 
1528  allocate (mesh%xmal1(mesh%num_nodes)); mesh%xmal1 = 0._f64
1529  allocate (mesh%xmal2(mesh%num_nodes)); mesh%xmal2 = 0._f64
1530  allocate (mesh%xmal3(mesh%num_nodes)); mesh%xmal3 = 0._f64
1531 
1532  allocate (mesh%nbcov(mesh%num_nodes + 1))
1533  allocate (mesh%nugcv(10*mesh%num_nodes))
1534  allocate (nuctfr(mesh%nmxfr))
1535 
1536  call poclis(mesh, ncotcu, nuctfr)
1537 
1538  deallocate (nuctfr)
1539 
1540 !quantities for boundary edges
1541 !
1542 ! Numbering of border dimensions in the direction
1543 ! trigonometric, by reference number, for
1544 ! carry out all diagnoses relating to
1545 ! borders.
1546 !
1547 ! kelfro - element number to which the dimension belongs
1548 ! kctfro - local number (1,2,3) of rating
1549 ! krefro - reference number of the dimension
1550 ! ksofro - numbers of the 2 vertices at the end of the side
1551 ! vnofro - normal vector components (inward)
1552 !
1553 ! nctfrt - Total number of border dimensions
1554 ! nctfro - Number of border dimensions by reference
1555 ! nctfrp - Array pointer (nb of dimension by ref)
1556 
1557 ! --- Initialisation (numerotation quelconque) -----------------
1558 
1559  allocate (mesh%nctfrp(0:mesh%nmxfr))
1560  allocate (mesh%nctfro(mesh%nmxfr))
1561  allocate (mesh%kctfro(mesh%nctfrt))
1562  allocate (mesh%kelfro(mesh%nctfrt))
1563  allocate (mesh%krefro(mesh%nctfrt))
1564  allocate (mesh%ksofro(2, mesh%nctfrt))
1565  allocate (mesh%vnofro(2, mesh%nctfrt))
1566 
1567  mesh%nctfro = 0
1568 
1569  ifr = 0
1570  do ict = 1, 3
1571  do iel = 1, mesh%num_triangles
1572  if (mesh%nvois(ict, iel) < 0) then
1573 
1574  iref = -mesh%nvois(ict, iel)
1575 
1576  ifr = ifr + 1
1577  mesh%kelfro(ifr) = iel !element
1578  mesh%kctfro(ifr) = ict !local edge
1579  mesh%krefro(ifr) = iref !reference edge
1580  mesh%ksofro(1, ifr) = mesh%nodes(ict, iel) ! 2 edge nodes
1581  mesh%ksofro(2, ifr) = mesh%nodes(mod(ict, 3) + 1, iel)
1582  mesh%nctfro(iref) = mesh%nctfro(iref) + 1 ! edge per reference
1583 
1584  end if
1585  end do
1586  end do
1587 
1588  mesh%nctfrp(0) = 0
1589  do ifr = 1, mesh%nmxfr
1590  mesh%nctfrp(ifr) = mesh%nctfrp(ifr - 1) + mesh%nctfro(ifr)
1591  end do
1592 
1593 !--- Loop over reference --------------------
1594 
1595  ictcl = 1
1596  do iref = 1, mesh%nmxfr
1597 
1598  nbcot = mesh%nctfro(iref)
1599  if (nbcot > 0) then
1600 
1601  ict1 = ictcl
1602 
1603  do ict = ict1, mesh%nctfrt
1604  jref = mesh%krefro(ict)
1605 
1606  if (jref == iref) then
1607 
1608  keltmp = mesh%kelfro(ict)
1609  kcttmp = mesh%kctfro(ict)
1610  kretmp = mesh%krefro(ict)
1611  ks1tmp = mesh%ksofro(1, ict)
1612  ks2tmp = mesh%ksofro(2, ict)
1613 
1614  mesh%kelfro(ict) = mesh%kelfro(ictcl)
1615  mesh%kctfro(ict) = mesh%kctfro(ictcl)
1616  mesh%krefro(ict) = mesh%krefro(ictcl)
1617  mesh%ksofro(1, ict) = mesh%ksofro(1, ictcl)
1618  mesh%ksofro(2, ict) = mesh%ksofro(2, ictcl)
1619 
1620  mesh%kelfro(ictcl) = keltmp
1621  mesh%kctfro(ictcl) = kcttmp
1622  mesh%krefro(ictcl) = kretmp
1623  mesh%ksofro(1, ictcl) = ks1tmp
1624  mesh%ksofro(2, ictcl) = ks2tmp
1625 
1626  ictcl = ictcl + 1
1627 
1628  end if
1629  end do
1630  end if
1631  end do
1632 
1633  do iref = 1, mesh%nmxfr
1634 
1635  nbcot = mesh%nctfro(iref)
1636 
1637  !One reference for the whole boundary
1638 
1639  if (nbcot == mesh%nctfrt) then
1640  ict1 = 1
1641 35 continue
1642  is2 = mesh%ksofro(2, ict1)
1643  do ict2 = ict1, mesh%nctfrt
1644  if (ict1 /= ict2) then
1645  is1 = mesh%ksofro(1, ict2)
1646  if (is1 == is2) then
1647 
1648  keltmp = mesh%kelfro(ict1 + 1)
1649  kcttmp = mesh%kctfro(ict1 + 1)
1650  kretmp = mesh%krefro(ict1 + 1)
1651  ks1tmp = mesh%ksofro(1, ict1 + 1)
1652  ks2tmp = mesh%ksofro(2, ict1 + 1)
1653 
1654  mesh%kelfro(ict1 + 1) = mesh%kelfro(ict2)
1655  mesh%kctfro(ict1 + 1) = mesh%kctfro(ict2)
1656  mesh%krefro(ict1 + 1) = mesh%krefro(ict2)
1657  mesh%ksofro(1, ict1 + 1) = mesh%ksofro(1, ict2)
1658  mesh%ksofro(2, ict1 + 1) = mesh%ksofro(2, ict2)
1659 
1660  mesh%kelfro(ict2) = keltmp
1661  mesh%kctfro(ict2) = kcttmp
1662  mesh%krefro(ict2) = kretmp
1663  mesh%ksofro(1, ict2) = ks1tmp
1664  mesh%ksofro(2, ict2) = ks2tmp
1665 
1666  if (ict1 < mesh%nctfrt) then
1667  ict1 = ict1 + 1
1668  goto 35
1669  end if
1670 
1671  end if
1672  end if
1673  end do
1674 
1675  ! ... several references for one boundary
1676 
1677  else if (nbcot > 1) then
1678 
1679  ictcl1 = mesh%nctfrp(iref - 1) + 1
1680  ictcl2 = mesh%nctfrp(iref)
1681 
1682  ict1 = ictcl1
1683 31 continue
1684  is1 = mesh%ksofro(1, ict1)
1685  do ict2 = ictcl1, ictcl2
1686  if (ict1 /= ict2) then
1687  is2 = mesh%ksofro(2, ict2)
1688  if (is1 == is2) then
1689  ict1 = ict2
1690  if (ict1 == ictcl1) then !Test si on tourne en rond ............
1691  goto 37
1692  end if
1693  goto 31
1694  end if
1695  end if
1696  end do
1697 
1698  keltmp = mesh%kelfro(ict1)
1699  kcttmp = mesh%kctfro(ict1)
1700  kretmp = mesh%krefro(ict1)
1701  ks1tmp = mesh%ksofro(1, ict1)
1702  ks2tmp = mesh%ksofro(2, ict1)
1703 
1704  mesh%kelfro(ict1) = mesh%kelfro(ictcl1)
1705  mesh%kctfro(ict1) = mesh%kctfro(ictcl1)
1706  mesh%krefro(ict1) = mesh%krefro(ictcl1)
1707 
1708  mesh%ksofro(1, ict1) = mesh%ksofro(1, ictcl1)
1709  mesh%ksofro(2, ict1) = mesh%ksofro(2, ictcl1)
1710 
1711  mesh%kelfro(ictcl1) = keltmp
1712  mesh%kctfro(ictcl1) = kcttmp
1713  mesh%krefro(ictcl1) = kretmp
1714 
1715  mesh%ksofro(1, ictcl1) = ks1tmp
1716  mesh%ksofro(2, ictcl1) = ks2tmp
1717 
1718 37 continue
1719  ict1 = ictcl1
1720 33 continue
1721  is2 = mesh%ksofro(2, ict1)
1722  do ict2 = ict1, ictcl2
1723  if (ict1 /= ict2) then
1724  is1 = mesh%ksofro(1, ict2)
1725  if (is1 == is2) then
1726 
1727  keltmp = mesh%kelfro(ict1 + 1)
1728  kcttmp = mesh%kctfro(ict1 + 1)
1729  kretmp = mesh%krefro(ict1 + 1)
1730  ks1tmp = mesh%ksofro(1, ict1 + 1)
1731  ks2tmp = mesh%ksofro(2, ict1 + 1)
1732 
1733  mesh%kelfro(ict1 + 1) = mesh%kelfro(ict2)
1734  mesh%kctfro(ict1 + 1) = mesh%kctfro(ict2)
1735  mesh%krefro(ict1 + 1) = mesh%krefro(ict2)
1736  mesh%ksofro(1, ict1 + 1) = mesh%ksofro(1, ict2)
1737  mesh%ksofro(2, ict1 + 1) = mesh%ksofro(2, ict2)
1738 
1739  mesh%kelfro(ict2) = keltmp
1740  mesh%kctfro(ict2) = kcttmp
1741  mesh%krefro(ict2) = kretmp
1742  mesh%ksofro(1, ict2) = ks1tmp
1743  mesh%ksofro(2, ict2) = ks2tmp
1744 
1745  ict1 = ict1 + 1
1746  goto 33
1747  end if
1748  end if
1749  end do
1750 
1751  if (ict1 < ictcl2) then
1752  ictcl1 = ict1 + 1
1753  ict1 = ictcl1
1754  goto 31
1755  end if
1756  end if
1757  end do
1758 
1759 !PN commented out because it is used only for particles simulation
1760 !! Modification de nvois ------------------------------------
1761 !! (Numero de reference change par le numero de cote)
1762 !
1763 !do ict=1,mesh%nctfrt
1764 ! ie=mesh%kelfro(ict)
1765 ! ic=mesh%kctfro(ict)
1766 ! mesh%nvois(ic,ie)=-ict
1767 !end do
1768 
1769 !--- Calcul des composantes du vecteur normal -----------------
1770 
1771  do ict = 1, mesh%nctfrt
1772 
1773  is1 = mesh%ksofro(1, ict)
1774  is2 = mesh%ksofro(2, ict)
1775 
1776  x1 = mesh%coord(1, is1)
1777  y1 = mesh%coord(2, is1)
1778  x2 = mesh%coord(1, is2)
1779  y2 = mesh%coord(2, is2)
1780 
1781  mesh%vnofro(1, ict) = -y2 + y1
1782  mesh%vnofro(2, ict) = x2 - x1
1783 
1784  end do
1785 
1786  mesh%analyzed = .true.
1787 
1788  end subroutine sll_s_analyze_triangular_mesh
1789 
1790 end module sll_m_triangular_meshes
sll_m_triangular_meshes::sll_o_create
Definition: sll_m_triangular_meshes.F90:116
sll_m_triangular_meshes::sll_o_delete
Definition: sll_m_triangular_meshes.F90:120
sll_m_triangular_meshes::sll_o_display
Definition: sll_m_triangular_meshes.F90:124
sll_m_triangular_meshes::sll_o_new_triangular_mesh_2d
Definition: sll_m_triangular_meshes.F90:128
sll_m_boundary_condition_descriptors
Describe different boundary conditions.
Definition: sll_m_boundary_condition_descriptors.F90:20
sll_m_constants
Fortran module where set some physical and mathematical constants.
Definition: sll_m_constants.F90:23
sll_m_constants::sll_p_pi
real(kind=f64), parameter, public sll_p_pi
Definition: sll_m_constants.F90:51
sll_m_hexagonal_meshes
Definition: sll_m_hexagonal_meshes.F90:23
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_triangular_meshes
Definition: sll_m_triangular_meshes.F90:19
sll_m_triangular_meshes::sll_s_write_triangular_mesh_mtv
subroutine, public sll_s_write_triangular_mesh_mtv(mesh, mtv_file)
Definition: sll_m_triangular_meshes.F90:604
sll_m_triangular_meshes::sll_s_analyze_triangular_mesh
subroutine, public sll_s_analyze_triangular_mesh(mesh)
Compute unstructured mesh quantities.
Definition: sll_m_triangular_meshes.F90:1466
sll_m_triangular_meshes::new_triangular_mesh_2d_from_square
type(sll_t_triangular_mesh_2d) function, pointer new_triangular_mesh_2d_from_square(nc_eta1, eta1_min, eta1_max, nc_eta2, eta2_min, eta2_max, bc)
Allocates the memory space for a new 2D triangular mesh on the heap,.
Definition: sll_m_triangular_meshes.F90:299
sll_m_triangular_meshes::new_triangular_mesh_2d_from_hex_mesh
type(sll_t_triangular_mesh_2d) function, pointer new_triangular_mesh_2d_from_hex_mesh(hex_mesh)
allocates the memory space for a new 2D triangular mesh on the heap, initializes it with the given he...
Definition: sll_m_triangular_meshes.F90:198
sll_m_triangular_meshes::display_triangular_mesh_2d
subroutine display_triangular_mesh_2d(mesh)
Displays mesh information on the terminal.
Definition: sll_m_triangular_meshes.F90:550
sll_m_triangular_meshes::sll_s_triangular_mesh_2d_init_from_file
subroutine, public sll_s_triangular_mesh_2d_init_from_file(m, maafil)
Initialize a new 2D triangular mesh.
Definition: sll_m_triangular_meshes.F90:183
sll_m_triangular_meshes::sll_s_read_from_file
subroutine, public sll_s_read_from_file(mesh, maafil)
Definition: sll_m_triangular_meshes.F90:873
sll_m_triangular_meshes::global_to_x1
real(kind=f64) function global_to_x1(mesh, i)
Definition: sll_m_triangular_meshes.F90:582
sll_m_triangular_meshes::initialize_triangular_mesh_2d
subroutine initialize_triangular_mesh_2d(mesh, nc_eta1, eta1_min, eta1_max, nc_eta2, eta2_min, eta2_max, bc)
Definition: sll_m_triangular_meshes.F90:374
sll_m_triangular_meshes::sll_s_triangular_mesh_2d_free
subroutine, public sll_s_triangular_mesh_2d_free(mesh)
Definition: sll_m_triangular_meshes.F90:865
sll_m_triangular_meshes::new_triangular_mesh_2d_from_file
type(sll_t_triangular_mesh_2d) function, pointer new_triangular_mesh_2d_from_file(maafil)
allocates the memory space for a new 2D triangular mesh on the heap, initializes it with the given ar...
Definition: sll_m_triangular_meshes.F90:167
sll_m_triangular_meshes::compute_areas
subroutine compute_areas(mesh)
Definition: sll_m_triangular_meshes.F90:1006
sll_m_triangular_meshes::get_cell_center
subroutine get_cell_center(mesh, iel, x1, x2)
Definition: sll_m_triangular_meshes.F90:942
sll_m_triangular_meshes::global_to_x2
real(kind=f64) function global_to_x2(mesh, i)
Definition: sll_m_triangular_meshes.F90:593
sll_m_triangular_meshes::sll_s_triangular_mesh_2d_init_from_square
subroutine, public sll_s_triangular_mesh_2d_init_from_square(m, nc_eta1, eta1_min, eta1_max, nc_eta2, eta2_min, eta2_max, bc)
Initialize a new 2D triangular mesh.
Definition: sll_m_triangular_meshes.F90:340
sll_m_triangular_meshes::poclis
subroutine poclis(mesh, ncotcu, nuctfr)
Definition: sll_m_triangular_meshes.F90:1277
sll_m_triangular_meshes::sll_s_map_to_circle
subroutine, public sll_s_map_to_circle(mesh, num_cells, order)
Map an hexagonal mesh on circle param[inout] mesh the triangular mesh built fron an hexagonal mesh pa...
Definition: sll_m_triangular_meshes.F90:973
sll_m_hexagonal_meshes::sll_t_hex_mesh_2d
2d hexagonal mesh
Definition: sll_m_hexagonal_meshes.F90:63
sll_m_triangular_meshes::sll_t_triangular_mesh_2d
2d hexagonal mesh
Definition: sll_m_triangular_meshes.F90:58
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1