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