Neko  0.9.99
A portable framework for high-order spectral element flow simulations
mesh.f90
Go to the documentation of this file.
1 ! Copyright (c) 2018-2023, The Neko Authors
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without
5 ! modification, are permitted provided that the following conditions
6 ! are met:
7 !
8 ! * Redistributions of source code must retain the above copyright
9 ! notice, this list of conditions and the following disclaimer.
10 !
11 ! * Redistributions in binary form must reproduce the above
12 ! copyright notice, this list of conditions and the following
13 ! disclaimer in the documentation and/or other materials provided
14 ! with the distribution.
15 !
16 ! * Neither the name of the authors nor the names of its
17 ! contributors may be used to endorse or promote products derived
18 ! from this software without specific prior written permission.
19 !
20 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 ! POSSIBILITY OF SUCH DAMAGE.
32 !
34 module mesh
35  use num_types, only : rp, dp, i8
36  use point, only : point_t
37  use element, only : element_t
38  use hex
39  use quad
40  use utils, only : neko_error, neko_warning
42  use tuple, only : tuple_i4_t, tuple4_i4_t
43  use htable
44  use datadist
45  use distdata
46  use comm
48  use math
49  use uset, only : uset_i8_t
50  use curve, only : curve_t
51  use logger, only : log_size
52  implicit none
53  private
54 
56  integer, public, parameter :: neko_msh_max_zlbls = 20
58  integer, public, parameter :: neko_msh_max_zlbl_len = 40
59 
60  type, private :: mesh_element_t
61  class(element_t), allocatable :: e
62  end type mesh_element_t
63 
64  type, public :: mesh_t
65  integer :: nelv
66  integer :: npts
67  integer :: gdim
68  integer :: mpts
69  integer :: mfcs
70  integer :: meds
71 
72  integer :: glb_nelv
73  integer :: glb_mpts
74  integer :: glb_mfcs
75  integer :: glb_meds
76 
77  integer :: offset_el
78  integer :: max_pts_id
79 
80  type(point_t), allocatable :: points(:)
81  type(mesh_element_t), allocatable :: elements(:)
82  logical, allocatable :: dfrmd_el(:)
83 
84  type(htable_i4_t) :: htp
85  type(htable_i4t4_t) :: htf
86  type(htable_i4t2_t) :: hte
87 
88  integer, allocatable :: facet_neigh(:,:)
89 
93  class(htable_t), allocatable :: facet_map
94  type(stack_i4_t), allocatable :: point_neigh(:)
95 
96  type(distdata_t) :: ddata
97  logical, allocatable :: neigh(:)
98  integer, allocatable :: neigh_order(:)
99 
100  integer(2), allocatable :: facet_type(:,:)
101 
103  type(facet_zone_t) :: inlet
104  type(facet_zone_t) :: outlet
105  type(facet_zone_t) :: outlet_normal
106  type(facet_zone_t) :: sympln
107  type(facet_zone_t), allocatable :: labeled_zones(:)
108  type(facet_zone_periodic_t) :: periodic
109  type(curve_t) :: curve
110 
111  logical :: lconn = .false.
112  logical :: ldist = .false.
113  logical :: lnumr = .false.
114  logical :: lgenc = .true.
115 
118  procedure(mesh_deform), pass(msh), pointer :: apply_deform => null()
119  contains
120  procedure, private, pass(this) :: init_nelv => mesh_init_nelv
121  procedure, private, pass(this) :: init_dist => mesh_init_dist
122  procedure, private, pass(this) :: add_quad => mesh_add_quad
123  procedure, private, pass(this) :: add_hex => mesh_add_hex
124  procedure, private, pass(this) :: add_edge => mesh_add_edge
125  procedure, private, pass(this) :: add_face => mesh_add_face
126  procedure, private, pass(this) :: add_point => mesh_add_point
127  procedure, private, pass(this) :: get_local_point => mesh_get_local_point
128  procedure, private, pass(this) :: get_local_edge => mesh_get_local_edge
129  procedure, private, pass(this) :: get_local_facet => mesh_get_local_facet
130  procedure, private, pass(this) :: get_global_edge => mesh_get_global_edge
131  procedure, private, pass(this) :: get_global_facet => mesh_get_global_facet
132  procedure, private, pass(this) :: is_shared_point => mesh_is_shared_point
133  procedure, private, pass(this) :: is_shared_edge => mesh_is_shared_edge
134  procedure, private, pass(this) :: is_shared_facet => mesh_is_shared_facet
135  procedure, pass(this) :: free => mesh_free
136  procedure, pass(this) :: finalize => mesh_finalize
137  procedure, pass(this) :: mark_wall_facet => mesh_mark_wall_facet
138  procedure, pass(this) :: mark_inlet_facet => mesh_mark_inlet_facet
139  procedure, pass(this) :: mark_outlet_facet => mesh_mark_outlet_facet
140  procedure, pass(this) :: mark_sympln_facet => mesh_mark_sympln_facet
141  procedure, pass(this) :: mark_periodic_facet => mesh_mark_periodic_facet
142  procedure, pass(this) :: mark_outlet_normal_facet => &
144  procedure, pass(this) :: mark_labeled_facet => mesh_mark_labeled_facet
145  procedure, pass(this) :: mark_curve_element => mesh_mark_curve_element
146  procedure, pass(this) :: apply_periodic_facet => mesh_apply_periodic_facet
147  procedure, pass(this) :: all_deformed => mesh_all_deformed
148  procedure, pass(this) :: get_facet_ids => mesh_get_facet_ids
149  procedure, pass(this) :: reset_periodic_ids => mesh_reset_periodic_ids
150  procedure, pass(this) :: create_periodic_ids => mesh_create_periodic_ids
151  procedure, pass(this) :: generate_conn => mesh_generate_conn
152  procedure, pass(this) :: have_point_glb_idx => mesh_have_point_glb_idx
154  generic :: init => init_nelv, init_dist
156  generic :: add_element => add_quad, add_hex
159  generic :: get_local => get_local_point, get_local_edge, get_local_facet
162  generic :: get_global => get_global_edge, get_global_facet
164  generic :: is_shared => is_shared_point, is_shared_edge, is_shared_facet
165  end type mesh_t
166 
167  abstract interface
168  subroutine mesh_deform(msh, x, y, z, lx, ly, lz)
169  import mesh_t
170  import rp
171  class(mesh_t) :: msh
172  integer, intent(in) :: lx, ly, lz
173  real(kind=rp), intent(inout) :: x(lx, ly, lz, msh%nelv)
174  real(kind=rp), intent(inout) :: y(lx, ly, lz, msh%nelv)
175  real(kind=rp), intent(inout) :: z(lx, ly, lz, msh%nelv)
176  end subroutine mesh_deform
177  end interface
178 
179  public :: mesh_deform
180 
181 contains
182 
184  subroutine mesh_init_nelv(this, gdim, nelv)
185  class(mesh_t), intent(inout) :: this
186  integer, intent(in) :: gdim
187  integer, intent(in) :: nelv
188  integer :: ierr
189  character(len=LOG_SIZE) :: log_buf
190 
191  call this%free()
192 
193  this%nelv = nelv
194  this%gdim = gdim
195 
196  if (this%nelv < 1) then
197  write(log_buf, '(A,I0,A)') 'MPI rank ', pe_rank, ' has zero elements'
198  call neko_warning(log_buf)
199  end if
200 
201  call mpi_allreduce(this%nelv, this%glb_nelv, 1, &
202  mpi_integer, mpi_sum, neko_comm, ierr)
203 
204  this%offset_el = 0
205  call mpi_exscan(this%nelv, this%offset_el, 1, &
206  mpi_integer, mpi_sum, neko_comm, ierr)
207 
208  call mesh_init_common(this)
209 
210  end subroutine mesh_init_nelv
211 
213  subroutine mesh_init_dist(this, gdim, dist)
214  class(mesh_t), intent(inout) :: this
215  integer, intent(in) :: gdim
216  type(linear_dist_t), intent(in) :: dist
217  character(len=LOG_SIZE) :: log_buf
218 
219  call this%free()
220 
221  this%nelv = dist%num_local()
222  if (this%nelv < 1) then
223  write(log_buf, '(A,I0,A)') 'MPI rank ', pe_rank, ' has zero elements'
224  call neko_warning(log_buf)
225  end if
226  this%glb_nelv = dist%num_global()
227  this%offset_el = dist%start_idx()
228  this%gdim = gdim
229 
230  call mesh_init_common(this)
231 
232  end subroutine mesh_init_dist
233 
234  subroutine mesh_init_common(this)
235  type(mesh_t), intent(inout) :: this
236  integer :: i
237  type(tuple_i4_t) :: facet_data
238 
239  this%max_pts_id = 0
240 
241  allocate(this%elements(this%nelv))
242  allocate(this%dfrmd_el(this%nelv))
243  if (this%gdim .eq. 3) then
244  do i = 1, this%nelv
245  allocate(hex_t::this%elements(i)%e)
246  end do
247  this%npts = neko_hex_npts
249  if (this%lgenc) then
250  allocate(htable_i4t4_t::this%facet_map)
251  select type (fmp => this%facet_map)
252  type is(htable_i4t4_t)
253  call fmp%init(this%nelv, facet_data)
254  end select
255 
256  allocate(this%facet_neigh(neko_hex_nfcs, this%nelv))
257 
258  call this%htf%init(this%nelv * neko_hex_nfcs, i)
259  call this%hte%init(this%nelv * neko_hex_neds, i)
260  end if
261  else if (this%gdim .eq. 2) then
262  do i = 1, this%nelv
263  allocate(quad_t::this%elements(i)%e)
264  end do
265  this%npts = neko_quad_npts
266  if (this%lgenc) then
267  allocate(htable_i4t2_t::this%facet_map)
268  select type (fmp => this%facet_map)
269  type is(htable_i4t2_t)
270  call fmp%init(this%nelv, facet_data)
271  end select
272 
273  allocate(this%facet_neigh(neko_quad_neds, this%nelv))
274 
275  call this%hte%init(this%nelv * neko_quad_neds, i)
276  end if
277  else
278  call neko_error("Invalid dimension")
279  end if
280 
282  allocate(this%points(this%npts*this%nelv))
283 
286  if (this%lgenc) then
287  allocate(this%point_neigh(this%gdim*this%npts*this%nelv))
288  do i = 1, this%gdim*this%npts*this%nelv
289  call this%point_neigh(i)%init()
290  end do
291  end if
292 
293  allocate(this%facet_type(2 * this%gdim, this%nelv))
294  this%facet_type = 0
295 
296  call this%htp%init(this%npts*this%nelv, i)
297 
298  call this%wall%init(this%nelv)
299  call this%inlet%init(this%nelv)
300  call this%outlet%init(this%nelv)
301  call this%outlet_normal%init(this%nelv)
302  call this%sympln%init(this%nelv)
303  call this%periodic%init(this%nelv)
304 
305  allocate(this%labeled_zones(neko_msh_max_zlbls))
306  do i = 1, neko_msh_max_zlbls
307  call this%labeled_zones(i)%init(this%nelv)
308  end do
309 
310  call this%curve%init(this%nelv)
311 
312  call distdata_init(this%ddata)
313 
314  allocate(this%neigh(0:pe_size-1))
315  this%neigh = .false.
316 
317  this%mpts = 0
318  this%mfcs = 0
319  this%meds = 0
320 
321  end subroutine mesh_init_common
322 
324  subroutine mesh_free(this)
325  class(mesh_t), intent(inout) :: this
326  integer :: i
327 
328  call this%htp%free()
329  call this%htf%free()
330  call this%hte%free()
331  call distdata_free(this%ddata)
332 
333 
334  if (allocated(this%dfrmd_el)) then
335  deallocate(this%dfrmd_el)
336  end if
337 
338  if (allocated(this%elements)) then
339  do i = 1, this%nelv
340  call this%elements(i)%e%free()
341  deallocate(this%elements(i)%e)
342  end do
343  deallocate(this%elements)
344  end if
345 
346  if (allocated(this%facet_map)) then
347  select type (fmp => this%facet_map)
348  type is(htable_i4t2_t)
349  call fmp%free()
350  type is(htable_i4t4_t)
351  call fmp%free()
352  end select
353  deallocate(this%facet_map)
354  end if
355 
356  if (allocated(this%facet_neigh)) then
357  deallocate(this%facet_neigh)
358  end if
359 
360  if (allocated(this%point_neigh)) then
361  do i = 1, this%gdim * this%npts * this%nelv
362  call this%point_neigh(i)%free()
363  end do
364  deallocate(this%point_neigh)
365  end if
366 
367  if (allocated(this%facet_type)) then
368  deallocate(this%facet_type)
369  end if
370  if (allocated(this%labeled_zones)) then
371  do i = 1, neko_msh_max_zlbls
372  call this%labeled_zones(i)%free()
373  end do
374  deallocate(this%labeled_zones)
375  end if
376 
377  if (allocated(this%neigh)) then
378  deallocate(this%neigh)
379  end if
380 
381  if (allocated(this%neigh_order)) then
382  deallocate(this%neigh_order)
383  end if
384 
385  if (allocated(this%points)) then
386  deallocate(this%points)
387  end if
388 
389  call this%wall%free()
390  call this%inlet%free()
391  call this%outlet%free()
392  call this%outlet_normal%free()
393  call this%sympln%free()
394  call this%periodic%free()
395  this%lconn = .false.
396  this%lnumr = .false.
397  this%ldist = .false.
398  this%lgenc = .true.
399 
400  end subroutine mesh_free
401 
402  subroutine mesh_finalize(this)
403  class(mesh_t), target, intent(inout) :: this
404  integer :: i
405 
406  call mesh_generate_flags(this)
407  call mesh_generate_conn(this)
408 
409  call this%wall%finalize()
410  call this%inlet%finalize()
411  call this%outlet%finalize()
412  call this%outlet_normal%finalize()
413  call this%sympln%finalize()
414  call this%periodic%finalize()
415  do i = 1, neko_msh_max_zlbls
416  call this%labeled_zones(i)%finalize()
417  end do
418  call this%curve%finalize()
419 
420  end subroutine mesh_finalize
421 
422  subroutine mesh_generate_flags(this)
423  type(mesh_t), intent(inout) :: this
424  real(kind=dp) :: u(3), v(3), w(3), temp
425  integer :: e
426 
427  do e = 1, this%nelv
428  if (this%gdim .eq. 2) then
429  this%dfrmd_el(e) = .false.
430  u = this%elements(e)%e%pts(2)%p%x - this%elements(e)%e%pts(1)%p%x
431  v = this%elements(e)%e%pts(3)%p%x - this%elements(e)%e%pts(1)%p%x
432  temp = u(1)*v(1) + u(2)*v(2)
433  if(.not. abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
434  else
435  this%dfrmd_el(e) = .false.
436  u = this%elements(e)%e%pts(2)%p%x - this%elements(e)%e%pts(1)%p%x
437  v = this%elements(e)%e%pts(3)%p%x - this%elements(e)%e%pts(1)%p%x
438  w = this%elements(e)%e%pts(5)%p%x - this%elements(e)%e%pts(1)%p%x
439  temp = u(1)*v(1) + u(2)*v(2) + u(3)*v(3)
440  if(.not. abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
441  temp = u(1)*w(1) + u(2)*w(2) + u(3)*w(3)
442  if(.not. abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
443  u = this%elements(e)%e%pts(7)%p%x - this%elements(e)%e%pts(8)%p%x
444  v = this%elements(e)%e%pts(6)%p%x - this%elements(e)%e%pts(8)%p%x
445  w = this%elements(e)%e%pts(4)%p%x - this%elements(e)%e%pts(8)%p%x
446  temp = u(1)*v(1) + u(2)*v(2) + u(3)*v(3)
447  if(.not. abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
448  temp = u(1)*w(1) + u(2)*w(2) + u(3)*w(3)
449  if(.not. abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
450  end if
451  end do
452  end subroutine mesh_generate_flags
453 
455  subroutine mesh_all_deformed(this)
456  class(mesh_t), intent(inout) :: this
457  this%dfrmd_el = .true.
458  end subroutine mesh_all_deformed
459 
461  subroutine mesh_generate_conn(this)
462  class(mesh_t), target, intent(inout) :: this
463  type(tuple_i4_t) :: edge
464  type(tuple4_i4_t) :: face, face_comp
465  type(tuple_i4_t) :: facet_data
466  type(stack_i4_t) :: neigh_order
467  class(element_t), pointer :: ep
468  type(tuple_i4_t) :: e
469  type(tuple4_i4_t) :: f
470  integer :: p_local_idx
471  integer :: el, id
472  integer :: i, j, k, ierr, el_glb_idx, n_sides, n_nodes, src, dst
473 
474  if (this%lconn) return
475 
476  if (.not. this%lgenc) return
477 
478  !If we generate connectivity, we do that here.
479  do el = 1, this%nelv
480  ep => this%elements(el)%e
481  select type(ep)
482  type is (hex_t)
483  do i = 1, neko_hex_npts
484  !Only for getting the id
485  call this%add_point(ep%pts(i)%p, id)
486  p_local_idx = this%get_local(this%points(id))
487  !should stack have inout on what we push? would be neat with in
488  id = ep%id()
489  call this%point_neigh(p_local_idx)%push(id)
490  end do
491  do i = 1, neko_hex_nfcs
492  call ep%facet_id(f, i)
493  call this%add_face(f)
494  end do
495 
496  do i = 1, neko_hex_neds
497  call ep%edge_id(e, i)
498  call this%add_edge(e)
499  end do
500  type is (quad_t)
501  do i = 1, neko_quad_npts
502  !Only for getting the id
503  call this%add_point(ep%pts(i)%p, id)
504  p_local_idx = this%get_local(this%points(id))
505  !should stack have inout on what we push? would be neat with in
506  id = ep%id()
507  call this%point_neigh(p_local_idx)%push(id)
508  end do
509 
510  do i = 1, neko_quad_neds
511  call ep%facet_id(e, i)
512  call this%add_edge(e)
513  end do
514  end select
515  end do
516 
517 
518  if (this%gdim .eq. 2) then
519  n_sides = 4
520  n_nodes = 2
521  else
522  n_sides = 6
523  n_nodes = 4
524  end if
525 
526  ! Compute global number of unique points
527  call mpi_allreduce(this%max_pts_id, this%glb_mpts, 1, &
528  mpi_integer, mpi_max, neko_comm, ierr)
529 
530  !
531  ! Find all (local) boundaries
532  !
533 
537  select type (fmp => this%facet_map)
538  type is(htable_i4t2_t)
539  do k = 1, 2
540  do i = 1, this%nelv
541  el_glb_idx = i + this%offset_el
542  do j = 1, n_sides
543  call this%elements(i)%e%facet_id(edge, j)
544 
545  ! Assume that all facets are on the exterior
546  facet_data%x = (/ 0, 0/)
547 
548  !check it this face has shown up earlier
549  if (fmp%get(edge, facet_data) .eq. 0) then
550  !if element is already recognized on face
551  if (facet_data%x(1) .eq. el_glb_idx ) then
552  this%facet_neigh(j, i) = facet_data%x(2)
553  else if( facet_data%x(2) .eq. el_glb_idx) then
554  this%facet_neigh(j, i) = facet_data%x(1)
555  !if this is the second element, arrange so low id is first
556  else if(facet_data%x(1) .gt. el_glb_idx) then
557  facet_data%x(2) = facet_data%x(1)
558  facet_data%x(1) = el_glb_idx
559  this%facet_neigh(j, i) = facet_data%x(2)
560  call fmp%set(edge, facet_data)
561  else if(facet_data%x(1) .lt. el_glb_idx) then
562  facet_data%x(2) = el_glb_idx
563  this%facet_neigh(j, i) = facet_data%x(1)
564  call fmp%set(edge, facet_data)
565  end if
566  else
567  facet_data%x(1) = el_glb_idx
568  this%facet_neigh(j, i) = facet_data%x(2)
569  call fmp%set(edge, facet_data)
570  end if
571  end do
572  end do
573  end do
574  type is(htable_i4t4_t)
575 
576  do k = 1, 2
577  do i = 1, this%nelv
578  el_glb_idx = i + this%offset_el
579  do j = 1, n_sides
580  call this%elements(i)%e%facet_id(face, j)
581 
582  facet_data%x = (/ 0, 0/)
583 
584  !check it this face has shown up earlier
585  if (fmp%get(face, facet_data) .eq. 0) then
586  !if element is already recognized on face
587  if (facet_data%x(1) .eq. el_glb_idx ) then
588  this%facet_neigh(j, i) = facet_data%x(2)
589  call this%elements(i)%e%facet_id(face_comp, &
590  j+(2*mod(j,2)-1))
591  if (face_comp .eq. face) then
592  facet_data%x(2) = el_glb_idx
593  this%facet_neigh(j, i) = facet_data%x(1)
594  call fmp%set(face, facet_data)
595  end if
596  else if( facet_data%x(2) .eq. el_glb_idx) then
597  this%facet_neigh(j, i) = facet_data%x(1)
598  !if this is the second element, arrange so low id is first
599  else if(facet_data%x(1) .gt. el_glb_idx) then
600  facet_data%x(2) = facet_data%x(1)
601  facet_data%x(1) = el_glb_idx
602  this%facet_neigh(j, i) = facet_data%x(2)
603  call fmp%set(face, facet_data)
604  else if(facet_data%x(1) .lt. el_glb_idx) then
605  facet_data%x(2) = el_glb_idx
606  this%facet_neigh(j, i) = facet_data%x(1)
607  call fmp%set(face, facet_data)
608  end if
609  else
610  facet_data%x(1) = el_glb_idx
611  this%facet_neigh(j, i) = 0
612  call fmp%set(face, facet_data)
613  end if
614  end do
615  end do
616  end do
617  class default
618  call neko_error('Invalid facet map')
619  end select
620 
621 
622  !
623  ! Find all external (between PEs) boundaries
624  !
625  if (pe_size .gt. 1) then
626 
628 
629  !
630  ! Generate neighbour exchange order
631  !
632  call neigh_order%init(pe_size)
633 
634  do i = 1, pe_size - 1
635  src = modulo(pe_rank - i + pe_size, pe_size)
636  dst = modulo(pe_rank + i, pe_size)
637  if (this%neigh(src) .or. this%neigh(dst)) then
638  j = i ! adhere to standards...
639  call neigh_order%push(j)
640  end if
641  end do
642 
643  allocate(this%neigh_order(neigh_order%size()))
644  select type(order => neigh_order%data)
645  type is (integer)
646  do i = 1, neigh_order%size()
647  this%neigh_order(i) = order(i)
648  end do
649  end select
650  call neigh_order%free()
651 
653  else
654  allocate(this%neigh_order(1))
655  this%neigh_order = 1
656  end if
657 
658  !
659  ! Find all internal/extenral edge connections
660  ! (Note it needs to be called after external point connections has
661  ! been established)
662  !
663  if (this%gdim .eq. 3) then
664  call mesh_generate_edge_conn(this)
665  end if
666 
667 
669 
670  this%lconn = .true.
671 
672  end subroutine mesh_generate_conn
673 
676  type(mesh_t), intent(inout) :: this
677  type(tuple_i4_t) :: edge, edge2
678  type(tuple4_i4_t) :: face, face2
679  type(tuple_i4_t) :: facet_data
680  type(stack_i4_t) :: buffer
681  type(mpi_status) :: status
682  type(mpi_request) :: send_req, recv_req
683  integer, allocatable :: recv_buffer(:)
684  integer :: i, j, k, el_glb_idx, n_sides, n_nodes, facet, element, l
685  integer :: max_recv, ierr, src, dst, n_recv, recv_side, neigh_el
686 
687 
688  if (this%gdim .eq. 2) then
689  n_sides = 4
690  n_nodes = 2
691  else
692  n_sides = 6
693  n_nodes = 4
694  end if
695 
696  call buffer%init()
697 
698  ! Build send buffers containing
699  ! [el_glb_idx, side number, facet_id (global ids of points)]
700  do i = 1, this%nelv
701  el_glb_idx = i + this%offset_el
702  do j = 1, n_sides
703  facet = j ! Adhere to standards...
704  if (this%facet_neigh(j, i) .eq. 0) then
705  if (n_nodes .eq. 2) then
706  call this%elements(i)%e%facet_id(edge, j)
707  call buffer%push(el_glb_idx)
708  call buffer%push(facet)
709  do k = 1, n_nodes
710  call buffer%push(edge%x(k))
711  end do
712  else
713  call this%elements(i)%e%facet_id(face, j)
714  call buffer%push(el_glb_idx)
715  call buffer%push(facet)
716  do k = 1, n_nodes
717  call buffer%push(face%x(k))
718  end do
719  end if
720  end if
721  end do
722  end do
723 
724 
725  call mpi_allreduce(buffer%size(), max_recv, 1, &
726  mpi_integer, mpi_max, neko_comm, ierr)
727 
728  allocate(recv_buffer(max_recv))
729 
730  do i = 1, size(this%neigh_order)
731  src = modulo(pe_rank - this%neigh_order(i) + pe_size, pe_size)
732  dst = modulo(pe_rank + this%neigh_order(i), pe_size)
733 
734  if (this%neigh(src)) then
735  call mpi_irecv(recv_buffer, max_recv, mpi_integer, &
736  src, 0, neko_comm, recv_req, ierr)
737  end if
738 
739  if (this%neigh(dst)) then
740  call mpi_isend(buffer%array(), buffer%size(), mpi_integer, &
741  dst, 0, neko_comm, send_req, ierr)
742  end if
743 
744  if (this%neigh(src)) then
745  call mpi_wait(recv_req, status, ierr)
746  call mpi_get_count(status, mpi_integer, n_recv, ierr)
747 
748  select type (fmp => this%facet_map)
749  type is(htable_i4t2_t)
750  do j = 1, n_recv, n_nodes + 2
751  neigh_el = recv_buffer(j)
752  recv_side = recv_buffer(j+1)
753 
754  edge = (/ recv_buffer(j+2), recv_buffer(j+3) /)
755 
756  facet_data = (/ 0, 0 /)
757  !Check if the face is present on this PE
758  if (fmp%get(edge, facet_data) .eq. 0) then
759  element = facet_data%x(1) - this%offset_el
760  !Check which side is connected
761  do l = 1, n_sides
762  call this%elements(element)%e%facet_id(edge2, l)
763  if(edge2 .eq. edge) then
764  facet = l
765  exit
766  end if
767  end do
768  this%facet_neigh(facet, element) = -neigh_el
769  facet_data%x(2) = -neigh_el
770 
771  ! Update facet map
772  call fmp%set(edge, facet_data)
773 
774  call distdata_set_shared_el_facet(this%ddata, element, facet)
775 
776  if (this%hte%get(edge, facet) .eq. 0) then
777  call distdata_set_shared_facet(this%ddata, facet)
778  else
779  call neko_error("Invalid shared edge")
780  end if
781 
782  end if
783 
784  end do
785  type is(htable_i4t4_t)
786  do j = 1, n_recv, n_nodes + 2
787  neigh_el = recv_buffer(j)
788  recv_side = recv_buffer(j+1)
789 
790  face%x = (/ recv_buffer(j+2), recv_buffer(j+3), &
791  recv_buffer(j+4), recv_buffer(j+5) /)
792 
793 
794  facet_data%x = (/ 0, 0 /)
795 
796  !Check if the face is present on this PE
797  if (fmp%get(face, facet_data) .eq. 0) then
798  ! Determine opposite side and update neighbor
799  element = facet_data%x(1) - this%offset_el
800  do l = 1, 6
801  call this%elements(element)%e%facet_id(face2, l)
802  if(face2 .eq. face) then
803  facet = l
804  exit
805  end if
806  end do
807  this%facet_neigh(facet, element) = -neigh_el
808  facet_data%x(2) = -neigh_el
809 
810  ! Update facet map
811  call fmp%set(face, facet_data)
812 
813  call distdata_set_shared_el_facet(this%ddata, element, facet)
814 
815  if (this%htf%get(face, facet) .eq. 0) then
816  call distdata_set_shared_facet(this%ddata, facet)
817  else
818  call neko_error("Invalid shared face")
819  end if
820 
821 
822  end if
823 
824  end do
825  end select
826  end if
827 
828  if (this%neigh(dst)) then
829  call mpi_wait(send_req, mpi_status_ignore, ierr)
830  end if
831 
832  end do
833 
834 
835  deallocate(recv_buffer)
836 
837  call buffer%free()
838 
839  end subroutine mesh_generate_external_facet_conn
840 
843  type(mesh_t), intent(inout) :: this
844  type(stack_i4_t) :: send_buffer
845  type(mpi_status) :: status
846  integer, allocatable :: recv_buffer(:)
847  integer :: i, j, k
848  integer :: max_recv, ierr, src, dst, n_recv, neigh_el
849  integer :: pt_glb_idx, pt_loc_idx, num_neigh
850  integer, contiguous, pointer :: neighs(:)
851 
852 
853  call send_buffer%init(this%mpts * 2)
854 
855  ! Build send buffers containing
856  ! [pt_glb_idx, #neigh, neigh id_1 ....neigh_id_n]
857  do i = 1, this%mpts
858  pt_glb_idx = this%points(i)%id() ! Adhere to standards...
859  num_neigh = this%point_neigh(i)%size()
860  call send_buffer%push(pt_glb_idx)
861  call send_buffer%push(num_neigh)
862 
863  neighs => this%point_neigh(i)%array()
864  do j = 1, this%point_neigh(i)%size()
865  call send_buffer%push(neighs(j))
866  end do
867  end do
868 
869  call mpi_allreduce(send_buffer%size(), max_recv, 1, &
870  mpi_integer, mpi_max, neko_comm, ierr)
871  allocate(recv_buffer(max_recv))
872 
873  do i = 1, pe_size - 1
874  src = modulo(pe_rank - i + pe_size, pe_size)
875  dst = modulo(pe_rank + i, pe_size)
876 
877  call mpi_sendrecv(send_buffer%array(), send_buffer%size(), &
878  mpi_integer, dst, 0, recv_buffer, max_recv, mpi_integer, src, 0, &
879  neko_comm, status, ierr)
880 
881  call mpi_get_count(status, mpi_integer, n_recv, ierr)
882 
883  j = 1
884  do while (j .le. n_recv)
885  pt_glb_idx = recv_buffer(j)
886  num_neigh = recv_buffer(j + 1)
887  ! Check if the point is present on this PE
888  pt_loc_idx = this%have_point_glb_idx(pt_glb_idx)
889  if (pt_loc_idx .gt. 0) then
890  this%neigh(src) = .true.
891  do k = 1, num_neigh
892  neigh_el = -recv_buffer(j + 1 + k)
893  call this%point_neigh(pt_loc_idx)%push(neigh_el)
894  call distdata_set_shared_point(this%ddata, pt_loc_idx)
895  end do
896  end if
897  j = j + (2 + num_neigh)
898  end do
899 
900  end do
901 
902  deallocate(recv_buffer)
903  call send_buffer%free()
904 
905  end subroutine mesh_generate_external_point_conn
906 
910  subroutine mesh_generate_edge_conn(this)
911  type(mesh_t), target, intent(inout) :: this
912  type(htable_iter_i4t2_t) :: it
913  type(tuple_i4_t), pointer :: edge
914  type(uset_i8_t), target :: edge_idx, ghost, owner
915  type(stack_i8_t), target :: send_buff
916  type(htable_i8_t) :: glb_to_loc
917  type(mpi_status) :: status
918  type(mpi_request) :: send_req, recv_req
919  integer, contiguous, pointer :: p1(:), p2(:), ns_id(:)
920  integer :: i, j, id, ierr, num_edge_glb, edge_offset, num_edge_loc
921  integer :: k, l , shared_offset, glb_nshared, n_glb_id
922  integer(kind=i8) :: C, glb_max, glb_id
923  integer(kind=i8), pointer :: glb_ptr
924  integer(kind=i8), allocatable :: recv_buff(:)
925  logical :: shared_edge
926  type(stack_i4_t), target :: non_shared_edges
927  integer :: max_recv, src, dst, n_recv
928 
929 
931  allocate(this%ddata%local_to_global_edge(this%meds))
932 
933  call edge_idx%init(this%hte%num_entries())
934  call send_buff%init(this%hte%num_entries())
935  call owner%init(this%hte%num_entries())
936 
937  call glb_to_loc%init(32, i)
938 
939  !
940  ! Determine/ constants used to generate unique global edge numbers
941  ! for shared edges
942  !
943  c = int(this%glb_nelv, i8) * int(neko_hex_neds, i8)
944 
945  num_edge_glb = 2* this%meds
946  call mpi_allreduce(mpi_in_place, num_edge_glb, 1, &
947  mpi_integer, mpi_sum, neko_comm, ierr)
948 
949  glb_max = int(num_edge_glb, i8)
950 
951  call non_shared_edges%init(this%hte%num_entries())
952 
953  call it%init(this%hte)
954  do while(it%next())
955  edge => it%key()
956  call it%data(id)
957 
958  k = this%have_point_glb_idx(edge%x(1))
959  l = this%have_point_glb_idx(edge%x(2))
960  p1 => this%point_neigh(k)%array()
961  p2 => this%point_neigh(l)%array()
962 
963  shared_edge = .false.
964 
965  ! Find edge neighbor from point neighbors
966  do i = 1, this%point_neigh(k)%size()
967  do j = 1, this%point_neigh(l)%size()
968  if ((p1(i) .eq. p2(j)) .and. &
969  (p1(i) .lt. 0) .and. (p2(j) .lt. 0)) then
970  call distdata_set_shared_edge(this%ddata, id)
971  shared_edge = .true.
972  end if
973  end do
974  end do
975 
976  ! Generate a unique id for the shared edge as,
977  ! ((e1 * C) + e2 )) + glb_max if e1 > e2
978  ! ((e2 * C) + e1 )) + glb_max if e2 > e1
979  if (shared_edge) then
980  glb_id = ((int(edge%x(1), i8)) + int(edge%x(2), i8)*c) + glb_max
981  call glb_to_loc%set(glb_id, id)
982  call edge_idx%add(glb_id)
983  call owner%add(glb_id) ! Always assume the PE is the owner
984  call send_buff%push(glb_id)
985  else
986  call non_shared_edges%push(id)
987  end if
988  end do
989 
990  ! Determine start offset for global numbering of locally owned edges
991  edge_offset = 0
992  num_edge_loc = non_shared_edges%size()
993  call mpi_exscan(num_edge_loc, edge_offset, 1, &
994  mpi_integer, mpi_sum, neko_comm, ierr)
995  edge_offset = edge_offset + 1
996 
997  ! Construct global numbering of locally owned edges
998  ns_id => non_shared_edges%array()
999  do i = 1, non_shared_edges%size()
1000  call distdata_set_local_to_global_edge(this%ddata, ns_id(i), edge_offset)
1001  edge_offset = edge_offset + 1
1002  end do
1003  nullify(ns_id)
1004 
1005  !
1006  ! Renumber shared edges into integer range
1007  !
1008 
1009  call mpi_allreduce(send_buff%size(), max_recv, 1, &
1010  mpi_integer, mpi_max, neko_comm, ierr)
1011 
1012  call ghost%init(send_buff%size())
1013 
1014  allocate(recv_buff(max_recv))
1015 
1016  do i = 1, size(this%neigh_order)
1017  src = modulo(pe_rank - this%neigh_order(i) + pe_size, pe_size)
1018  dst = modulo(pe_rank + this%neigh_order(i), pe_size)
1019 
1020  if (this%neigh(src)) then
1021  call mpi_irecv(recv_buff, max_recv, mpi_integer8, &
1022  src, 0, neko_comm, recv_req, ierr)
1023  end if
1024 
1025  if (this%neigh(dst)) then
1026  ! We should use the %array() procedure, which works great for
1027  ! GNU, Intel and NEC, but it breaks horribly on Cray when using
1028  ! certain data types
1029  select type(sbarray=>send_buff%data)
1030  type is (integer(i8))
1031  call mpi_isend(sbarray, send_buff%size(), mpi_integer8, &
1032  dst, 0, neko_comm, send_req, ierr)
1033  end select
1034  end if
1035 
1036  if (this%neigh(src)) then
1037  call mpi_wait(recv_req, status, ierr)
1038  call mpi_get_count(status, mpi_integer8, n_recv, ierr)
1039 
1040  do j = 1, n_recv
1041  if ((edge_idx%element(recv_buff(j))) .and. (src .lt. pe_rank)) then
1042  call ghost%add(recv_buff(j))
1043  call owner%remove(recv_buff(j))
1044  end if
1045  end do
1046  end if
1047 
1048  if (this%neigh(dst)) then
1049  call mpi_wait(send_req, mpi_status_ignore, ierr)
1050  end if
1051  end do
1052 
1053 
1054  ! Determine start offset for global numbering of shared edges
1055  glb_nshared = num_edge_loc
1056  call mpi_allreduce(mpi_in_place, glb_nshared, 1, &
1057  mpi_integer, mpi_sum, neko_comm, ierr)
1058 
1059  shared_offset = 0
1060  call mpi_exscan(owner%size(), shared_offset, 1, &
1061  mpi_integer, mpi_sum, neko_comm, ierr)
1062  shared_offset = shared_offset + glb_nshared + 1
1063 
1064  ! Renumber locally owned set of shared edges
1065  call send_buff%clear()
1066  call owner%iter_init()
1067  do while (owner%iter_next())
1068  glb_ptr => owner%iter_value()
1069  if (glb_to_loc%get(glb_ptr, id) .eq. 0) then
1070  call distdata_set_local_to_global_edge(this%ddata, id, shared_offset)
1071 
1072  ! Add new number to send data as [old_glb_id new_glb_id] for each edge
1073  call send_buff%push(glb_ptr) ! Old glb_id integer*8
1074  glb_id = int(shared_offset, i8) ! Waste some space here...
1075  call send_buff%push(glb_id) ! New glb_id integer*4
1076 
1077  shared_offset = shared_offset + 1
1078  else
1079  call neko_error('Invalid edge id')
1080  end if
1081  end do
1082  nullify(glb_ptr)
1083 
1084  ! Determine total number of unique edges in the mesh
1085  ! (This can probably be done in a clever way...)
1086  this%glb_meds = shared_offset -1
1087  call mpi_allreduce(mpi_in_place, this%glb_meds, 1, &
1088  mpi_integer, mpi_max, neko_comm, ierr)
1089 
1090  !
1091  ! Update ghosted edges with new global id
1092  !
1093 
1094  call mpi_allreduce(send_buff%size(), max_recv, 1, &
1095  mpi_integer, mpi_max, neko_comm, ierr)
1096 
1097  deallocate(recv_buff)
1098  allocate(recv_buff(max_recv))
1099 
1100 
1101  do i = 1, size(this%neigh_order)
1102  src = modulo(pe_rank - this%neigh_order(i) + pe_size, pe_size)
1103  dst = modulo(pe_rank + this%neigh_order(i), pe_size)
1104 
1105  if (this%neigh(src)) then
1106  call mpi_irecv(recv_buff, max_recv, mpi_integer8, &
1107  src, 0, neko_comm, recv_req, ierr)
1108  end if
1109 
1110  if (this%neigh(dst)) then
1111  ! We should use the %array() procedure, which works great for
1112  ! GNU, Intel and NEC, but it breaks horribly on Cray when using
1113  ! certain data types
1114  select type(sbarray=>send_buff%data)
1115  type is (integer(i8))
1116  call mpi_isend(sbarray, send_buff%size(), mpi_integer8, &
1117  dst, 0, neko_comm, send_req, ierr)
1118  end select
1119  end if
1120 
1121  if (this%neigh(src)) then
1122  call mpi_wait(recv_req, status, ierr)
1123  call mpi_get_count(status, mpi_integer8, n_recv, ierr)
1124 
1125  do j = 1, n_recv, 2
1126  if (ghost%element(recv_buff(j))) then
1127  if (glb_to_loc%get(recv_buff(j), id) .eq. 0) then
1128  n_glb_id = int(recv_buff(j + 1 ), 4)
1129  call distdata_set_local_to_global_edge(this%ddata, id, &
1130  n_glb_id)
1131  else
1132  call neko_error('Invalid edge id')
1133  end if
1134  end if
1135  end do
1136  end if
1137 
1138  if (this%neigh(dst)) then
1139  call mpi_wait(send_req, mpi_status_ignore, ierr)
1140  end if
1141  end do
1142 
1143  deallocate(recv_buff)
1144  call glb_to_loc%free()
1145  call send_buff%free()
1146  call edge_idx%free()
1147  call non_shared_edges%free()
1148  call ghost%free()
1149  call owner%free()
1150 
1151  end subroutine mesh_generate_edge_conn
1152 
1155  type(mesh_t), target, intent(inout) :: this
1156  type(htable_iter_i4t4_t), target :: face_it
1157  type(htable_iter_i4t2_t), target :: edge_it
1158  type(tuple4_i4_t), pointer :: face, fd(:)
1159  type(tuple_i4_t), pointer :: edge, ed(:)
1160  type(tuple_i4_t) :: facet_data
1161  type(tuple4_i4_t) :: recv_face
1162  type(tuple_i4_t) :: recv_edge
1163  type(stack_i4t4_t) :: face_owner
1164  type(htable_i4t4_t) :: face_ghost
1165  type(stack_i4t2_t) :: edge_owner
1166  type(htable_i4t2_t) :: edge_ghost
1167  type(stack_i4_t) :: send_buff
1168  type(mpi_status) :: status
1169  type(mpi_request) :: send_req, recv_req
1170  integer, allocatable :: recv_buff(:)
1171  integer :: non_shared_facets, shared_facets, facet_offset
1172  integer :: id, glb_nshared, shared_offset, owned_facets
1173  integer :: i, j, ierr, max_recv, src, dst, n_recv
1174 
1175  shared_facets = this%ddata%shared_facet%size()
1176 
1178  if (this%gdim .eq. 2) then
1179  allocate(this%ddata%local_to_global_facet(this%meds))
1180  call edge_owner%init(this%meds)
1181  call edge_ghost%init(64, i)
1182  non_shared_facets = this%hte%num_entries() - shared_facets
1183  else
1184  allocate(this%ddata%local_to_global_facet(this%mfcs))
1185  call face_owner%init(this%mfcs)
1186  call face_ghost%init(64, i)
1187  non_shared_facets = this%htf%num_entries() - shared_facets
1188  end if
1189 
1191 
1192  facet_offset = 0
1193  call mpi_exscan(non_shared_facets, facet_offset, 1, &
1194  mpi_integer, mpi_sum, neko_comm, ierr)
1195  facet_offset = facet_offset + 1
1196 
1197  ! Determine ownership of shared facets
1198  if (this%gdim .eq. 2) then
1199  call edge_it%init(this%hte)
1200  do while (edge_it%next())
1201  call edge_it%data(id)
1202  edge => edge_it%key()
1203  if (.not. this%ddata%shared_facet%element(id)) then
1204  call distdata_set_local_to_global_facet(this%ddata, &
1205  id, facet_offset)
1206  facet_offset = facet_offset + 1
1207  else
1208  select type(fmp => this%facet_map)
1209  type is(htable_i4t2_t)
1210  if (fmp%get(edge, facet_data) .eq. 0) then
1211  if (facet_data%x(2) .lt. 0) then
1212  if (abs(facet_data%x(2)) .lt. (this%offset_el + 1)) then
1213  call edge_ghost%set(edge, id)
1214  else
1215  call edge_owner%push(edge)
1216  end if
1217  else
1218  call neko_error("Invalid edge neigh.")
1219  end if
1220  end if
1221  end select
1222  end if
1223  end do
1224  owned_facets = edge_owner%size()
1225  else
1226  call face_it%init(this%htf)
1227  do while (face_it%next())
1228  call face_it%data(id)
1229  face => face_it%key()
1230  if (.not. this%ddata%shared_facet%element(id)) then
1231  call distdata_set_local_to_global_facet(this%ddata, &
1232  id, facet_offset)
1233  facet_offset = facet_offset + 1
1234  else
1235  select type(fmp => this%facet_map)
1236  type is(htable_i4t4_t)
1237  if (fmp%get(face, facet_data) .eq. 0) then
1238  if (facet_data%x(2) .lt. 0) then
1239  if (abs(facet_data%x(2)) .lt. (this%offset_el + 1)) then
1240  call face_ghost%set(face, id)
1241  else
1242  call face_owner%push(face)
1243  end if
1244  else
1245  call neko_error("Invalid face neigh.")
1246  end if
1247  end if
1248  end select
1249  end if
1250  end do
1251  owned_facets = face_owner%size()
1252  end if
1253 
1254  ! Determine start offset for global numbering of shared facets
1255  glb_nshared = non_shared_facets
1256  call mpi_allreduce(mpi_in_place, glb_nshared, 1, &
1257  mpi_integer, mpi_sum, neko_comm, ierr)
1258 
1259  shared_offset = 0
1260  call mpi_exscan(owned_facets, shared_offset, 1, &
1261  mpi_integer, mpi_sum, neko_comm, ierr)
1262  shared_offset = shared_offset + glb_nshared + 1
1263 
1264  if (this%gdim .eq. 2) then
1265 
1266  if (owned_facets .gt. 32) then
1267  call send_buff%init(owned_facets)
1268  else
1269  call send_buff%init()
1270  end if
1271 
1272  ed => edge_owner%array()
1273  do i = 1, edge_owner%size()
1274  if (this%hte%get(ed(i), id) .eq. 0) then
1275  call distdata_set_local_to_global_facet(this%ddata, id, &
1276  shared_offset)
1277 
1278  ! Add new number to send buffer
1279  ! [edge id1 ... edge idn new_glb_id]
1280  do j = 1, 2
1281  call send_buff%push(ed(i)%x(j))
1282  end do
1283  call send_buff%push(shared_offset)
1284 
1285  shared_offset = shared_offset + 1
1286  end if
1287  end do
1288 
1289  else
1290 
1291  if (owned_facets .gt. 32) then
1292  call send_buff%init(owned_facets)
1293  else
1294  call send_buff%init()
1295  end if
1296 
1297  fd => face_owner%array()
1298  do i = 1, face_owner%size()
1299  if (this%htf%get(fd(i), id) .eq. 0) then
1300  call distdata_set_local_to_global_facet(this%ddata, id, &
1301  shared_offset)
1302 
1303  ! Add new number to send buffer
1304  ! [face id1 ... face idn new_glb_id]
1305  do j = 1, 4
1306  call send_buff%push(fd(i)%x(j))
1307  end do
1308  call send_buff%push(shared_offset)
1309 
1310  shared_offset = shared_offset + 1
1311  end if
1312  end do
1313  nullify(fd)
1314 
1315  end if
1316 
1317  ! Determine total number of unique facets in the mesh
1318  ! (This can probably be done in a clever way...)
1319  this%glb_mfcs = shared_offset - 1
1320  call mpi_allreduce(mpi_in_place, this%glb_mfcs, 1, &
1321  mpi_integer, mpi_max, neko_comm, ierr)
1322 
1323  !
1324  ! Update ghosted facets with new global id
1325  !
1326 
1327  call mpi_allreduce(send_buff%size(), max_recv, 1, &
1328  mpi_integer, mpi_max, neko_comm, ierr)
1329 
1330  allocate(recv_buff(max_recv))
1331 
1333  do i = 1, size(this%neigh_order)
1334  src = modulo(pe_rank - this%neigh_order(i) + pe_size, pe_size)
1335  dst = modulo(pe_rank + this%neigh_order(i), pe_size)
1336 
1337  if (this%neigh(src)) then
1338  call mpi_irecv(recv_buff, max_recv, mpi_integer, &
1339  src, 0, neko_comm, recv_req, ierr)
1340  end if
1341 
1342  if (this%neigh(dst)) then
1343  call mpi_isend(send_buff%array(), send_buff%size(), mpi_integer, &
1344  dst, 0, neko_comm, send_req, ierr)
1345  end if
1346 
1347  if (this%neigh(src)) then
1348  call mpi_wait(recv_req, status, ierr)
1349  call mpi_get_count(status, mpi_integer, n_recv, ierr)
1350 
1351  if (this%gdim .eq. 2) then
1352  do j = 1, n_recv, 3
1353 
1354  recv_edge = (/recv_buff(j), recv_buff(j+1)/)
1355 
1356  ! Check if the PE has the shared edge
1357  if (edge_ghost%get(recv_edge, id) .eq. 0) then
1358  call distdata_set_local_to_global_facet(this%ddata, &
1359  id, recv_buff(j+2))
1360  end if
1361  end do
1362  else
1363  do j = 1, n_recv, 5
1364 
1365  recv_face = (/recv_buff(j), recv_buff(j+1), &
1366  recv_buff(j+2), recv_buff(j+3) /)
1367 
1368  ! Check if the PE has the shared face
1369  if (face_ghost%get(recv_face, id) .eq. 0) then
1370  call distdata_set_local_to_global_facet(this%ddata, &
1371  id, recv_buff(j+4))
1372  end if
1373  end do
1374  end if
1375  end if
1376 
1377  if (this%neigh(dst)) then
1378  call mpi_wait(send_req, mpi_status_ignore, ierr)
1379  end if
1380 
1381  end do
1382 
1383  if (this%gdim .eq. 2) then
1384  call edge_owner%free()
1385  call edge_ghost%free()
1386  else
1387  call face_owner%free()
1388  call face_ghost%free()
1389  end if
1390 
1391  call send_buff%free()
1392  deallocate(recv_buff)
1393 
1394  end subroutine mesh_generate_facet_numbering
1395 
1396 
1398  subroutine mesh_add_quad(this, el, p1, p2, p3, p4)
1399  class(mesh_t), target, intent(inout) :: this
1400  integer, value :: el
1401  type(point_t), target, intent(inout) :: p1, p2, p3, p4
1402  class(element_t), pointer :: ep
1403  integer :: p(4), el_glb_idx
1404  type(tuple_i4_t) :: e
1405 
1406  ! Connectivity invalidated if a new element is added
1407  this%lconn = .false.
1408 
1409  ! Numbering invalidated if a new element is added
1410  this%lnumr = .false.
1411 
1412  call this%add_point(p1, p(1))
1413  call this%add_point(p2, p(2))
1414  call this%add_point(p3, p(3))
1415  call this%add_point(p4, p(4))
1416 
1417  ep => this%elements(el)%e
1418  el_glb_idx = el + this%offset_el
1419 
1420  select type(ep)
1421  type is (quad_t)
1422  call ep%init(el_glb_idx, &
1423  this%points(p(1)), this%points(p(2)), &
1424  this%points(p(3)), this%points(p(4)))
1425 
1426 
1427  class default
1428  call neko_error('Invalid element type')
1429  end select
1430 
1431  end subroutine mesh_add_quad
1432 
1434  subroutine mesh_add_hex(this, el, p1, p2, p3, p4, p5, p6, p7, p8)
1435  class(mesh_t), target, intent(inout) :: this
1436  integer, value :: el
1437  type(point_t), target, intent(inout) :: p1, p2, p3, p4, p5, p6, p7, p8
1438  class(element_t), pointer :: ep
1439  integer :: p(8), el_glb_idx
1440  type(tuple4_i4_t) :: f
1441  type(tuple_i4_t) :: e
1442 
1443  ! Connectivity invalidated if a new element is added
1444  this%lconn = .false.
1445 
1446  ! Numbering invalidated if a new element is added
1447  this%lnumr = .false.
1448 
1449  call this%add_point(p1, p(1))
1450  call this%add_point(p2, p(2))
1451  call this%add_point(p3, p(3))
1452  call this%add_point(p4, p(4))
1453  call this%add_point(p5, p(5))
1454  call this%add_point(p6, p(6))
1455  call this%add_point(p7, p(7))
1456  call this%add_point(p8, p(8))
1457 
1458  ep => this%elements(el)%e
1459  el_glb_idx = el + this%offset_el
1460  select type(ep)
1461  type is (hex_t)
1462  call ep%init(el_glb_idx, &
1463  this%points(p(1)), this%points(p(2)), &
1464  this%points(p(3)), this%points(p(4)), &
1465  this%points(p(5)), this%points(p(6)), &
1466  this%points(p(7)), this%points(p(8)))
1467  class default
1468  call neko_error('Invalid element type')
1469  end select
1470 
1471  end subroutine mesh_add_hex
1472 
1474  subroutine mesh_add_point(this, p, idx)
1475  class(mesh_t), intent(inout) :: this
1476  type(point_t), intent(inout) :: p
1477  integer, intent(inout) :: idx
1478  integer :: tmp
1479 
1480  tmp = p%id()
1481 
1482  this%max_pts_id = max(this%max_pts_id, tmp)
1483 
1484  if (tmp .le. 0) then
1485  call neko_error("Invalid point id")
1486  end if
1487 
1488  if (this%htp%get(tmp, idx) .gt. 0) then
1489  this%mpts = this%mpts + 1
1490  call this%htp%set(tmp, this%mpts)
1491  this%points(this%mpts) = p
1492  idx = this%mpts
1493  end if
1494 
1495  end subroutine mesh_add_point
1496 
1498  subroutine mesh_add_face(this, f)
1499  class(mesh_t), intent(inout) :: this
1500  type(tuple4_i4_t), intent(inout) :: f
1501  integer :: idx
1502 
1503  if (this%htf%get(f, idx) .gt. 0) then
1504  this%mfcs = this%mfcs + 1
1505  call this%htf%set(f, this%mfcs)
1506  end if
1507 
1508  end subroutine mesh_add_face
1509 
1511  subroutine mesh_add_edge(this, e)
1512  class(mesh_t), intent(inout) :: this
1513  type(tuple_i4_t), intent(inout) :: e
1514  integer :: idx
1515 
1516  if (this%hte%get(e, idx) .gt. 0) then
1517  this%meds = this%meds + 1
1518  call this%hte%set(e, this%meds)
1519  end if
1520 
1521  end subroutine mesh_add_edge
1522 
1524  subroutine mesh_mark_wall_facet(this, f, e)
1525  class(mesh_t), intent(inout) :: this
1526  integer, intent(inout) :: f
1527  integer, intent(inout) :: e
1528 
1529  if (e .gt. this%nelv) then
1530  call neko_error('Invalid element index')
1531  end if
1532 
1533  if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1534  (this%gdim .eq. 3 .and. f .gt. 6)) then
1535  call neko_error('Invalid facet index')
1536  end if
1537  this%facet_type(f, e) = 2
1538  call this%wall%add_facet(f, e)
1539 
1540  end subroutine mesh_mark_wall_facet
1541 
1543  subroutine mesh_mark_curve_element(this, e, curve_data, curve_type)
1544  class(mesh_t), intent(inout) :: this
1545  integer, intent(in) :: e
1546  real(kind=dp), dimension(5,12), intent(in) :: curve_data
1547  integer, dimension(12), intent(in) :: curve_type
1548 
1549  if (e .gt. this%nelv) then
1550  call neko_error('Invalid element index')
1551  end if
1552  if ((this%gdim .eq. 2 .and. sum(curve_type(5:8)) .gt. 0) ) then
1553  call neko_error('Invalid curve element')
1554  end if
1555  call this%curve%add_element(e, curve_data, curve_type)
1556 
1557  end subroutine mesh_mark_curve_element
1558 
1559 
1561  subroutine mesh_mark_inlet_facet(this, f, e)
1562  class(mesh_t), intent(inout) :: this
1563  integer, intent(in) :: f
1564  integer, intent(in) :: e
1565 
1566  if (e .gt. this%nelv) then
1567  call neko_error('Invalid element index')
1568  end if
1569 
1570  if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1571  (this%gdim .eq. 3 .and. f .gt. 6)) then
1572  call neko_error('Invalid facet index')
1573  end if
1574  this%facet_type(f, e) = 2
1575  call this%inlet%add_facet(f, e)
1576 
1577  end subroutine mesh_mark_inlet_facet
1578 
1580  subroutine mesh_mark_labeled_facet(this, f, e, label)
1581  class(mesh_t), intent(inout) :: this
1582  integer, intent(in) :: f
1583  integer, intent(in) :: e
1584  integer, intent(in) :: label
1585 
1586  if (e .gt. this%nelv) then
1587  call neko_error('Invalid element index')
1588  end if
1589 
1590  if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1591  (this%gdim .eq. 3 .and. f .gt. 6)) then
1592  call neko_error('Invalid facet index')
1593  end if
1594  call this%labeled_zones(label)%add_facet(f, e)
1595  this%facet_type(f,e) = -label
1596 
1597  end subroutine mesh_mark_labeled_facet
1598 
1599 
1601  subroutine mesh_mark_outlet_normal_facet(this, f, e)
1602  class(mesh_t), intent(inout) :: this
1603  integer, intent(inout) :: f
1604  integer, intent(inout) :: e
1605 
1606  if (e .gt. this%nelv) then
1607  call neko_error('Invalid element index')
1608  end if
1609 
1610  if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1611  (this%gdim .eq. 3 .and. f .gt. 6)) then
1612  call neko_error('Invalid facet index')
1613  end if
1614  this%facet_type(f, e) = 1
1615  call this%outlet_normal%add_facet(f, e)
1616 
1617  end subroutine mesh_mark_outlet_normal_facet
1618 
1619 
1621  subroutine mesh_mark_outlet_facet(this, f, e)
1622  class(mesh_t), intent(inout) :: this
1623  integer, intent(inout) :: f
1624  integer, intent(inout) :: e
1625 
1626  if (e .gt. this%nelv) then
1627  call neko_error('Invalid element index')
1628  end if
1629 
1630  if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1631  (this%gdim .eq. 3 .and. f .gt. 6)) then
1632  call neko_error('Invalid facet index')
1633  end if
1634  this%facet_type(f, e) = 1
1635  call this%outlet%add_facet(f, e)
1636 
1637  end subroutine mesh_mark_outlet_facet
1638 
1640  subroutine mesh_mark_sympln_facet(this, f, e)
1641  class(mesh_t), intent(inout) :: this
1642  integer, intent(inout) :: f
1643  integer, intent(inout) :: e
1644 
1645  if (e .gt. this%nelv) then
1646  call neko_error('Invalid element index')
1647  end if
1648 
1649  if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1650  (this%gdim .eq. 3 .and. f .gt. 6)) then
1651  call neko_error('Invalid facet index')
1652  end if
1653  this%facet_type(f, e) = 2
1654  call this%sympln%add_facet(f, e)
1655 
1656  end subroutine mesh_mark_sympln_facet
1657 
1659  subroutine mesh_mark_periodic_facet(this, f, e, pf, pe, pids)
1660  class(mesh_t), intent(inout) :: this
1661  integer, intent(in) :: f
1662  integer, intent(in) :: e
1663  integer, intent(in) :: pf
1664  integer, intent(in) :: pe
1665  integer, intent(inout) :: pids(4)
1666  integer, dimension(4) :: org_ids
1667 
1668  call this%get_facet_ids(f, e, org_ids)
1669  call this%periodic%add_periodic_facet(f, e, pf, pe, pids, org_ids)
1670  end subroutine mesh_mark_periodic_facet
1671 
1673  subroutine mesh_get_facet_ids(this, f, e, pids)
1674  class(mesh_t), intent(inout) :: this
1675  integer, intent(in) :: f
1676  integer, intent(in) :: e
1677  integer, intent(inout) :: pids(4)
1678  type(point_t), pointer :: pi
1679  type(tuple4_i4_t) :: t
1680  type(tuple_i4_t) :: t2
1681 
1682  select type(ele => this%elements(e)%e)
1683  type is(hex_t)
1684  call ele%facet_order(t,f)
1685  pids = t%x
1686  type is(quad_t)
1687  call ele%facet_order(t2,f)
1688  pids(1) = t2%x(1)
1689  pids(2) = t2%x(2)
1690  pids(3) = 0
1691  pids(4) = 0
1692  end select
1693  end subroutine mesh_get_facet_ids
1694 
1696  subroutine mesh_reset_periodic_ids(this)
1697  class(mesh_t), intent(inout) :: this
1698  integer :: i,j
1699  integer :: f
1700  integer :: e
1701  integer :: pf
1702  integer :: pe
1703  integer :: org_ids(4), pids(4)
1704  type(point_t), pointer :: pi
1705  integer, dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1706  2,6,8,4,&
1707  1,2,6,5,&
1708  3,4,8,7,&
1709  1,2,4,3,&
1710  5,6,8,7/),&
1711  (/4,6/))
1712  integer, dimension(2, 4) :: edge_nodes = reshape((/1,3,&
1713  2,4,&
1714  1,2,&
1715  3,4 /),&
1716  (/2,4/))
1717 
1718  do i = 1, this%periodic%size
1719  e = this%periodic%facet_el(i)%x(2)
1720  f = this%periodic%facet_el(i)%x(1)
1721  pe = this%periodic%p_facet_el(i)%x(2)
1722  pf = this%periodic%p_facet_el(i)%x(1)
1723  pids = this%periodic%p_ids(i)%x
1724  call this%get_facet_ids(f, e, pids)
1725  this%periodic%p_ids(i)%x = pids
1726  end do
1727  do i = 1, this%periodic%size
1728  e = this%periodic%facet_el(i)%x(2)
1729  f = this%periodic%facet_el(i)%x(1)
1730  org_ids = this%periodic%org_ids(i)%x
1731  select type(ele => this%elements(e)%e)
1732  type is(hex_t)
1733  do j = 1, 4
1734  pi => ele%pts(face_nodes(j,f))%p
1735  call pi%set_id(org_ids(j))
1736  end do
1737  type is(quad_t)
1738  do j = 1, 2
1739  pi => ele%pts(edge_nodes(j,f))%p
1740  call pi%set_id(org_ids(j))
1741  end do
1742  end select
1743  end do
1744  end subroutine mesh_reset_periodic_ids
1745 
1747  subroutine mesh_create_periodic_ids(this, f, e, pf, pe)
1748  class(mesh_t), intent(inout) :: this
1749  integer, intent(in) :: f
1750  integer, intent(in) :: e
1751  integer, intent(in) :: pf
1752  integer, intent(in) :: pe
1753  type(point_t), pointer :: pi, pj
1754  real(kind=dp) :: l(3)
1755  integer :: i, j, id, p_local_idx, match
1756  type(tuple4_i4_t) :: ft
1757  type(tuple_i4_t) :: et
1758  integer, dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1759  2,6,8,4,&
1760  1,2,6,5,&
1761  3,4,8,7,&
1762  1,2,4,3,&
1763  5,6,8,7/),&
1764  (/4,6/))
1765  integer, dimension(2, 4) :: edge_nodes = reshape((/1,3,&
1766  2,4,&
1767  1,2,&
1768  3,4 /),&
1769  (/2,4/))
1770 
1771  select type(ele => this%elements(e)%e)
1772  type is(hex_t)
1773  select type(elp => this%elements(pe)%e)
1774  type is(hex_t)
1775  l = 0d0
1776  do i = 1, 4
1777  l = l + ele%pts(face_nodes(i,f))%p%x(1:3) - &
1778  elp%pts(face_nodes(i,pf))%p%x(1:3)
1779  end do
1780  l = l/4
1781  do i = 1, 4
1782  pi => ele%pts(face_nodes(i,f))%p
1783  match = 0
1784  do j = 1, 4
1785  pj => elp%pts(face_nodes(j,pf))%p
1786  if (norm2(pi%x(1:3) - pj%x(1:3) - l) .lt. 1d-7) then
1787  id = min(pi%id(), pj%id())
1788  call pi%set_id(id)
1789  call pj%set_id(id)
1790  p_local_idx = this%get_local(this%points(id))
1791  match = match + 1
1792  end if
1793  end do
1794  if ( match .gt. 1) then
1795  call neko_error('Multiple matches when creating periodic ids')
1796  else if (match .eq. 0) then
1797  call neko_error('Cannot find matching periodic point')
1798  end if
1799  end do
1800  end select
1801  type is(quad_t)
1802  select type(elp => this%elements(pe)%e)
1803  type is(quad_t)
1804  l = 0d0
1805  do i = 1, 2
1806  l = l + ele%pts(edge_nodes(i,f))%p%x(1:3) - &
1807  elp%pts(edge_nodes(i,pf))%p%x(1:3)
1808  end do
1809  l = l/2
1810  do i = 1, 2
1811  pi => ele%pts(edge_nodes(i,f))%p
1812  do j = 1, 2
1813  pj => elp%pts(edge_nodes(j,pf))%p
1814  !whatabout thie tolerance?
1815  if (norm2(pi%x(1:3) - pj%x(1:3) - l) .lt. 1d-7) then
1816  id = min(pi%id(), pj%id())
1817  call pi%set_id(id)
1818  call pj%set_id(id)
1819  p_local_idx = this%get_local(this%points(id))
1820  end if
1821  end do
1822  end do
1823  end select
1824  end select
1825  end subroutine mesh_create_periodic_ids
1826 
1829  subroutine mesh_apply_periodic_facet(this, f, e, pf, pe, pids)
1830  class(mesh_t), intent(inout) :: this
1831  integer, intent(in) :: f
1832  integer, intent(in) :: e
1833  integer, intent(in) :: pf
1834  integer, intent(in) :: pe
1835  integer, intent(inout) :: pids(4)
1836  type(point_t), pointer :: pi
1837  integer :: i, id, p_local_idx
1838  type(tuple4_i4_t) :: ft
1839  type(tuple_i4_t) :: et
1840  integer, dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1841  2,6,8,4,&
1842  1,2,6,5,&
1843  3,4,8,7,&
1844  1,2,4,3,&
1845  5,6,8,7/),&
1846  (/4,6/))
1847  select type(ele => this%elements(e)%e)
1848  type is(hex_t)
1849  do i = 1, 4
1850  pi => ele%pts(face_nodes(i,f))%p
1851  call pi%set_id(pids(i))
1852  call this%add_point(pi, id)
1853  p_local_idx = this%get_local(this%points(id))
1854  end do
1855  end select
1856 
1857  end subroutine mesh_apply_periodic_facet
1858 
1860  function mesh_get_local_point(this, p) result(local_id)
1861  class(mesh_t), intent(inout) :: this
1862  type(point_t), intent(inout) :: p
1863  integer :: local_id
1864  integer :: tmp
1865 
1867  tmp = p%id()
1868 
1869  if (this%htp%get(tmp, local_id) .gt. 0) then
1870  call neko_error('Invalid global id (local point)')
1871  end if
1872 
1873  end function mesh_get_local_point
1874 
1877  function mesh_get_local_edge(this, e) result(local_id)
1878  class(mesh_t), intent(inout) :: this
1879  type(tuple_i4_t), intent(inout) :: e
1880  integer :: local_id
1881 
1882  if (this%hte%get(e, local_id) .gt. 0) then
1883  call neko_error('Invalid global id (local edge)')
1884  end if
1885 
1886  end function mesh_get_local_edge
1887 
1889  function mesh_get_local_facet(this, f) result(local_id)
1890  class(mesh_t), intent(inout) :: this
1891  type(tuple4_i4_t), intent(inout) :: f
1892  integer :: local_id
1893 
1894  if (this%htf%get(f, local_id) .gt. 0) then
1895  call neko_error('Invalid global id (local facet)')
1896  end if
1897 
1898  end function mesh_get_local_facet
1899 
1901  function mesh_get_global_edge(this, e) result(global_id)
1902  class(mesh_t), intent(inout) :: this
1903  type(tuple_i4_t), intent(inout) :: e
1904  integer :: global_id
1905 
1906  global_id = this%get_local(e)
1907 
1908  if (this%gdim .eq. 2) then
1909  if (pe_size .gt. 1) then
1910  global_id = this%ddata%local_to_global_facet(global_id)
1911  end if
1912  else
1913  if (pe_size .gt. 1) then
1914  global_id = this%ddata%local_to_global_edge(global_id)
1915  end if
1916  end if
1917 
1918  end function mesh_get_global_edge
1919 
1921  function mesh_get_global_facet(this, f) result(global_id)
1922  class(mesh_t), intent(inout) :: this
1923  type(tuple4_i4_t), intent(inout) :: f
1924  integer :: global_id
1925 
1926  global_id = this%get_local_facet(f)
1927 
1928  if (pe_size .gt. 1) then
1929  global_id = this%ddata%local_to_global_facet(global_id)
1930  end if
1931 
1932  end function mesh_get_global_facet
1933 
1934 
1938  function mesh_have_point_glb_idx(this, index) result(local_id)
1939  class(mesh_t), intent(inout) :: this
1940  integer, intent(inout) :: index
1941  integer :: local_id
1942 
1943  if (this%htp%get(index, local_id) .eq. 1) then
1944  local_id = -1
1945  end if
1946 
1947  end function mesh_have_point_glb_idx
1948 
1949 
1951  function mesh_is_shared_point(this, p) result(shared)
1952  class(mesh_t), intent(inout) :: this
1953  type(point_t), intent(inout) :: p
1954  integer :: local_index
1955  logical shared
1956 
1957  local_index = this%get_local(p)
1958  shared = this%ddata%shared_point%element(local_index)
1959 
1960  end function mesh_is_shared_point
1961 
1962 
1965  function mesh_is_shared_edge(this, e) result(shared)
1966  class(mesh_t), intent(inout) :: this
1967  type(tuple_i4_t), intent(inout) :: e
1968  integer :: local_index
1969  logical shared
1970  local_index = this%get_local(e)
1971  if (this%gdim .eq. 2) then
1972  shared = this%ddata%shared_facet%element(local_index)
1973  else
1974  shared = this%ddata%shared_edge%element(local_index)
1975  end if
1976  end function mesh_is_shared_edge
1977 
1979  function mesh_is_shared_facet(this, f) result(shared)
1980  class(mesh_t), intent(inout) :: this
1981  type(tuple4_i4_t), intent(inout) :: f
1982  integer :: local_index
1983  logical shared
1984 
1985  local_index = this%get_local(f)
1986  shared = this%ddata%shared_facet%element(local_index)
1987 
1988  end function mesh_is_shared_facet
1989 
1990 end module mesh
Definition: comm.F90:1
integer pe_rank
MPI rank.
Definition: comm.F90:28
type(mpi_comm) neko_comm
MPI communicator.
Definition: comm.F90:16
integer pe_size
MPI size of communicator.
Definition: comm.F90:31
Defines a domain as a subset of facets in a mesh.
Definition: curve.f90:2
Defines practical data distributions.
Definition: datadist.f90:34
Distributed mesh data.
Definition: distdata.f90:34
subroutine, public distdata_set_shared_point(ddata, point)
Mark a point as shared.
Definition: distdata.f90:127
subroutine, public distdata_set_shared_el_facet(ddata, element, side)
Mark an element's facet as shared.
Definition: distdata.f90:96
subroutine, public distdata_set_local_to_global_facet(ddata, local, global)
Set local to global mapping (facets)
Definition: distdata.f90:136
subroutine, public distdata_free(ddata)
Free a distdata type.
Definition: distdata.f90:76
subroutine, public distdata_set_local_to_global_edge(ddata, local, global)
Set local to global mapping (edges)
Definition: distdata.f90:146
subroutine, public distdata_init(ddata)
Initialise a distdata type.
Definition: distdata.f90:62
subroutine, public distdata_set_shared_facet(ddata, facet)
Mark a facet as shared.
Definition: distdata.f90:108
subroutine, public distdata_set_shared_edge(ddata, edge)
Mark an element's edge as shared.
Definition: distdata.f90:118
Defines a zone as a subset of facets in a mesh.
Definition: facet_zone.f90:34
Defines a hexahedron element.
Definition: hex.f90:34
integer, parameter, public neko_hex_npts
Number of points.
Definition: hex.f90:42
integer, parameter, public neko_hex_nfcs
Number of faces.
Definition: hex.f90:43
integer, parameter, public neko_hex_neds
Number of edges.
Definition: hex.f90:44
Implements a hash table ADT.
Definition: htable.f90:36
Logging routines.
Definition: log.f90:34
integer, parameter, public log_size
Definition: log.f90:42
Definition: math.f90:60
Defines a mesh.
Definition: mesh.f90:34
subroutine mesh_generate_flags(this)
Definition: mesh.f90:423
subroutine mesh_generate_facet_numbering(this)
Generate a unique facet numbering.
Definition: mesh.f90:1155
subroutine mesh_init_dist(this, gdim, dist)
Initialise a mesh this based on a distribution dist.
Definition: mesh.f90:214
logical function mesh_is_shared_point(this, p)
Check if a point is shared.
Definition: mesh.f90:1952
subroutine mesh_mark_wall_facet(this, f, e)
Mark facet f in element e as a wall.
Definition: mesh.f90:1525
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
Definition: mesh.f90:56
subroutine mesh_mark_sympln_facet(this, f, e)
Mark facet f in element e as a symmetry plane.
Definition: mesh.f90:1641
subroutine mesh_add_face(this, f)
Add a unique face represented as a 4-tuple to the mesh.
Definition: mesh.f90:1499
integer function mesh_get_global_edge(this, e)
Return the global id of an edge e.
Definition: mesh.f90:1902
subroutine mesh_generate_edge_conn(this)
Generate element-element connectivity via edges both between internal and between PEs.
Definition: mesh.f90:911
subroutine mesh_add_hex(this, el, p1, p2, p3, p4, p5, p6, p7, p8)
Add a hexahedral element to the mesh this.
Definition: mesh.f90:1435
subroutine mesh_add_edge(this, e)
Add a unique edge represented as a 2-tuple to the mesh.
Definition: mesh.f90:1512
subroutine mesh_add_point(this, p, idx)
Add a unique point to the mesh.
Definition: mesh.f90:1475
subroutine mesh_free(this)
Deallocate a mesh this.
Definition: mesh.f90:325
subroutine mesh_mark_labeled_facet(this, f, e, label)
Mark facet f in element e with label.
Definition: mesh.f90:1581
integer, parameter, public neko_msh_max_zlbl_len
Max length of a zone label.
Definition: mesh.f90:58
subroutine mesh_mark_periodic_facet(this, f, e, pf, pe, pids)
Mark facet f in element e as periodic with (pf, pe)
Definition: mesh.f90:1660
integer function mesh_get_local_edge(this, e)
Return the local id of an edge e.
Definition: mesh.f90:1878
subroutine mesh_mark_outlet_normal_facet(this, f, e)
Mark facet f in element e as an outlet normal.
Definition: mesh.f90:1602
subroutine mesh_all_deformed(this)
Set all elements as if they are deformed.
Definition: mesh.f90:456
logical function mesh_is_shared_edge(this, e)
Check if an edge is shared.
Definition: mesh.f90:1966
subroutine mesh_mark_inlet_facet(this, f, e)
Mark facet f in element e as an inlet.
Definition: mesh.f90:1562
subroutine mesh_init_common(this)
Definition: mesh.f90:235
subroutine mesh_get_facet_ids(this, f, e, pids)
Get original ids of periodic points.
Definition: mesh.f90:1674
subroutine mesh_create_periodic_ids(this, f, e, pf, pe)
Creates common ids for matching periodic points.
Definition: mesh.f90:1748
subroutine mesh_reset_periodic_ids(this)
Reset ids of periodic points to their original ids.
Definition: mesh.f90:1697
subroutine mesh_init_nelv(this, gdim, nelv)
Initialise a mesh this with nelv elements.
Definition: mesh.f90:185
logical function mesh_is_shared_facet(this, f)
Check if a facet is shared.
Definition: mesh.f90:1980
integer function mesh_get_global_facet(this, f)
Return the local id of a face f.
Definition: mesh.f90:1922
integer function mesh_have_point_glb_idx(this, index)
Check if the mesh has a point given its global index.
Definition: mesh.f90:1939
subroutine mesh_generate_external_point_conn(this)
Generate element-element connectivity via points between PEs.
Definition: mesh.f90:843
subroutine mesh_apply_periodic_facet(this, f, e, pf, pe, pids)
Replaces the periodic point's id with a common id for matching periodic points.
Definition: mesh.f90:1830
subroutine mesh_finalize(this)
Definition: mesh.f90:403
subroutine mesh_generate_external_facet_conn(this)
Generate element-element connectivity via facets between PEs.
Definition: mesh.f90:676
subroutine mesh_mark_curve_element(this, e, curve_data, curve_type)
Mark element e as a curve element.
Definition: mesh.f90:1544
subroutine mesh_generate_conn(this)
Generate element-to-element connectivity.
Definition: mesh.f90:462
subroutine mesh_mark_outlet_facet(this, f, e)
Mark facet f in element e as an outlet.
Definition: mesh.f90:1622
integer function mesh_get_local_facet(this, f)
Return the local id of a face f.
Definition: mesh.f90:1890
integer function mesh_get_local_point(this, p)
Return the local id of a point p.
Definition: mesh.f90:1861
subroutine mesh_add_quad(this, el, p1, p2, p3, p4)
Add a quadrilateral element to the mesh this.
Definition: mesh.f90:1399
integer, parameter, public i8
Definition: num_types.f90:7
integer, parameter, public dp
Definition: num_types.f90:9
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Implements a point.
Definition: point.f90:35
Defines a quadrilateral element.
Definition: quad.f90:34
integer, parameter, public neko_quad_neds
Number of edges.
Definition: quad.f90:43
integer, parameter, public neko_quad_npts
Number of points.
Definition: quad.f90:42
Implements a dynamic stack ADT.
Definition: stack.f90:35
Implements a n-tuple.
Definition: tuple.f90:34
Implements an unordered set ADT.
Definition: uset.f90:35
Utilities.
Definition: utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition: utils.f90:245
Defines wall boundary conditions.
Definition: wall.f90:34
Load-balanced linear distribution .
Definition: datadist.f90:50
Base type for an element.
Definition: element.f90:44
Hexahedron element.
Definition: hex.f90:63
Integer based hash table.
Definition: htable.f90:82
Integer 2-tuple based hash table.
Definition: htable.f90:122
Integer 4-tuple based hash table.
Definition: htable.f90:132
Integer*8 based hash table.
Definition: htable.f90:92
Iterator for an integer based 2-tuple hash table.
Definition: htable.f90:202
Iterator for an integer based 4-tuple hash table.
Definition: htable.f90:211
Base type for a hash table.
Definition: htable.f90:55
A point in with coordinates .
Definition: point.f90:43
Quadrilateral element.
Definition: quad.f90:58
Integer based stack.
Definition: stack.f90:63
Integer 2-tuple based stack.
Definition: stack.f90:84
Integer 4-tuple based stack.
Definition: stack.f90:91
Integer*8 based stack.
Definition: stack.f90:70
Integer based 4-tuple.
Definition: tuple.f90:69
Integer based 2-tuple.
Definition: tuple.f90:51
Integer*8 based unordered set.
Definition: uset.f90:74
#define max(a, b)
Definition: tensor.cu:40