81 logical,
allocatable :: dfrmd_el(:)
87 integer,
allocatable :: facet_neigh(:,:)
96 logical,
allocatable :: neigh(:)
97 integer,
allocatable :: neigh_order(:)
99 integer(2),
allocatable :: facet_type(:,:)
110 logical :: lconn = .false.
111 logical :: ldist = .false.
112 logical :: lnumr = .false.
113 logical :: lgenc = .true.
117 procedure(
mesh_deform), pass(msh),
pointer :: apply_deform => null()
141 procedure, pass(this) :: mark_outlet_normal_facet => &
153 generic :: init => init_nelv, init_dist
155 generic :: add_element => add_quad, add_hex
158 generic :: get_local => get_local_point, get_local_edge, get_local_facet
161 generic :: get_global => get_global_edge, get_global_facet
163 generic :: is_shared => is_shared_point, is_shared_edge, is_shared_facet
171 integer,
intent(in) :: lx, ly, lz
172 real(kind=
rp),
intent(inout) :: x(lx, ly, lz, msh%nelv)
173 real(kind=
rp),
intent(inout) :: y(lx, ly, lz, msh%nelv)
174 real(kind=
rp),
intent(inout) :: z(lx, ly, lz, msh%nelv)
184 class(
mesh_t),
intent(inout) :: this
185 integer,
intent(in) :: gdim
186 integer,
intent(in) :: nelv
194 call mpi_allreduce(this%nelv, this%glb_nelv, 1, &
198 call mpi_exscan(this%nelv, this%offset_el, 1, &
207 class(
mesh_t),
intent(inout) :: this
208 integer,
intent(in) :: gdim
213 this%nelv = dist%num_local()
214 this%glb_nelv = dist%num_global()
215 this%offset_el = dist%start_idx()
223 type(
mesh_t),
intent(inout) :: this
229 allocate(this%elements(this%nelv))
230 allocate(this%dfrmd_el(this%nelv))
231 if (this%gdim .eq. 3)
then
233 allocate(
hex_t::this%elements(i)%e)
239 select type (fmp => this%facet_map)
241 call fmp%init(this%nelv, facet_data)
249 else if (this%gdim .eq. 2)
then
251 allocate(
quad_t::this%elements(i)%e)
256 select type (fmp => this%facet_map)
258 call fmp%init(this%nelv, facet_data)
270 allocate(this%points(this%npts*this%nelv))
275 allocate(this%point_neigh(this%gdim*this%npts*this%nelv))
276 do i = 1, this%gdim*this%npts*this%nelv
277 call this%point_neigh(i)%init()
281 allocate(this%facet_type(2 * this%gdim, this%nelv))
284 call this%htp%init(this%npts*this%nelv, i)
286 call this%wall%init(this%nelv)
287 call this%inlet%init(this%nelv)
288 call this%outlet%init(this%nelv)
289 call this%outlet_normal%init(this%nelv)
290 call this%sympln%init(this%nelv)
291 call this%periodic%init(this%nelv)
295 call this%labeled_zones(i)%init(this%nelv)
298 call this%curve%init(this%nelv)
302 allocate(this%neigh(0:
pe_size-1))
313 class(
mesh_t),
intent(inout) :: this
322 if (
allocated(this%dfrmd_el))
then
323 deallocate(this%dfrmd_el)
326 if (
allocated(this%elements))
then
328 call this%elements(i)%e%free()
329 deallocate(this%elements(i)%e)
331 deallocate(this%elements)
334 if (
allocated(this%facet_map))
then
335 select type (fmp => this%facet_map)
341 deallocate(this%facet_map)
344 if (
allocated(this%facet_neigh))
then
345 deallocate(this%facet_neigh)
348 if (
allocated(this%point_neigh))
then
349 do i = 1, this%gdim * this%npts * this%nelv
350 call this%point_neigh(i)%free()
352 deallocate(this%point_neigh)
355 if (
allocated(this%facet_type))
then
356 deallocate(this%facet_type)
358 if (
allocated(this%labeled_zones))
then
360 call this%labeled_zones(i)%free()
362 deallocate(this%labeled_zones)
365 if (
allocated(this%neigh))
then
366 deallocate(this%neigh)
369 if (
allocated(this%neigh_order))
then
370 deallocate(this%neigh_order)
373 if (
allocated(this%points))
then
374 deallocate(this%points)
377 call this%wall%free()
378 call this%inlet%free()
379 call this%outlet%free()
380 call this%outlet_normal%free()
381 call this%sympln%free()
382 call this%periodic%free()
391 class(
mesh_t),
target,
intent(inout) :: this
397 call this%wall%finalize()
398 call this%inlet%finalize()
399 call this%outlet%finalize()
400 call this%outlet_normal%finalize()
401 call this%sympln%finalize()
402 call this%periodic%finalize()
404 call this%labeled_zones(i)%finalize()
406 call this%curve%finalize()
411 type(
mesh_t),
intent(inout) :: this
412 real(kind=
dp) :: u(3), v(3), w(3), temp
416 if (this%gdim .eq. 2)
then
417 this%dfrmd_el(e) = .false.
418 u = this%elements(e)%e%pts(2)%p%x - this%elements(e)%e%pts(1)%p%x
419 v = this%elements(e)%e%pts(3)%p%x - this%elements(e)%e%pts(1)%p%x
420 temp = u(1)*v(1) + u(2)*v(2)
421 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
423 this%dfrmd_el(e) = .false.
424 u = this%elements(e)%e%pts(2)%p%x - this%elements(e)%e%pts(1)%p%x
425 v = this%elements(e)%e%pts(3)%p%x - this%elements(e)%e%pts(1)%p%x
426 w = this%elements(e)%e%pts(5)%p%x - this%elements(e)%e%pts(1)%p%x
427 temp = u(1)*v(1) + u(2)*v(2) + u(3)*v(3)
428 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
429 temp = u(1)*w(1) + u(2)*w(2) + u(3)*w(3)
430 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
431 u = this%elements(e)%e%pts(7)%p%x - this%elements(e)%e%pts(8)%p%x
432 v = this%elements(e)%e%pts(6)%p%x - this%elements(e)%e%pts(8)%p%x
433 w = this%elements(e)%e%pts(4)%p%x - this%elements(e)%e%pts(8)%p%x
434 temp = u(1)*v(1) + u(2)*v(2) + u(3)*v(3)
435 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
436 temp = u(1)*w(1) + u(2)*w(2) + u(3)*w(3)
437 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
444 class(
mesh_t),
intent(inout) :: this
445 this%dfrmd_el = .true.
450 class(
mesh_t),
target,
intent(inout) :: this
458 integer :: p_local_idx, res
460 integer :: i, j, k, ierr, el_glb_idx, n_sides, n_nodes, src, dst
462 if (this%lconn)
return
464 if (.not. this%lgenc)
return
468 ep => this%elements(el)%e
473 call this%add_point(ep%pts(i)%p, id)
474 p_local_idx = this%get_local(this%points(id))
477 call this%point_neigh(p_local_idx)%push(id)
480 call ep%facet_id(f, i)
481 call this%add_face(f)
485 call ep%edge_id(e, i)
486 call this%add_edge(e)
491 call this%add_point(ep%pts(i)%p, id)
492 p_local_idx = this%get_local(this%points(id))
495 call this%point_neigh(p_local_idx)%push(id)
499 call ep%facet_id(e, i)
500 call this%add_edge(e)
506 if (this%gdim .eq. 2)
then
515 call mpi_allreduce(this%max_pts_id, this%glb_mpts, 1, &
525 select type (fmp => this%facet_map)
529 el_glb_idx = i + this%offset_el
531 call this%elements(i)%e%facet_id(edge, j)
534 facet_data%x = (/ 0, 0/)
537 if (fmp%get(edge, facet_data) .eq. 0)
then
539 if (facet_data%x(1) .eq. el_glb_idx )
then
540 this%facet_neigh(j, i) = facet_data%x(2)
541 else if( facet_data%x(2) .eq. el_glb_idx)
then
542 this%facet_neigh(j, i) = facet_data%x(1)
544 else if(facet_data%x(1) .gt. el_glb_idx)
then
545 facet_data%x(2) = facet_data%x(1)
546 facet_data%x(1) = el_glb_idx
547 this%facet_neigh(j, i) = facet_data%x(2)
548 call fmp%set(edge, facet_data)
549 else if(facet_data%x(1) .lt. el_glb_idx)
then
550 facet_data%x(2) = el_glb_idx
551 this%facet_neigh(j, i) = facet_data%x(1)
552 call fmp%set(edge, facet_data)
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)
566 el_glb_idx = i + this%offset_el
568 call this%elements(i)%e%facet_id(face, j)
570 facet_data%x = (/ 0, 0/)
573 if (fmp%get(face, facet_data) .eq. 0)
then
575 if (facet_data%x(1) .eq. el_glb_idx )
then
576 this%facet_neigh(j, i) = facet_data%x(2)
577 call this%elements(i)%e%facet_id(face_comp, &
579 if (face_comp .eq. face)
then
580 facet_data%x(2) = el_glb_idx
581 this%facet_neigh(j, i) = facet_data%x(1)
582 call fmp%set(face, facet_data)
584 else if( facet_data%x(2) .eq. el_glb_idx)
then
585 this%facet_neigh(j, i) = facet_data%x(1)
587 else if(facet_data%x(1) .gt. el_glb_idx)
then
588 facet_data%x(2) = facet_data%x(1)
589 facet_data%x(1) = el_glb_idx
590 this%facet_neigh(j, i) = facet_data%x(2)
591 call fmp%set(face, facet_data)
592 else if(facet_data%x(1) .lt. el_glb_idx)
then
593 facet_data%x(2) = el_glb_idx
594 this%facet_neigh(j, i) = facet_data%x(1)
595 call fmp%set(face, facet_data)
598 facet_data%x(1) = el_glb_idx
599 this%facet_neigh(j, i) = 0
600 call fmp%set(face, facet_data)
625 if (this%neigh(src) .or. this%neigh(dst))
then
627 call neigh_order%push(j)
631 allocate(this%neigh_order(neigh_order%size()))
632 select type(order => neigh_order%data)
634 do i = 1, neigh_order%size()
635 this%neigh_order(i) = order(i)
638 call neigh_order%free()
642 allocate(this%neigh_order(1))
651 if (this%gdim .eq. 3)
then
664 type(
mesh_t),
intent(inout) :: this
669 type(mpi_status) :: status
670 type(mpi_request) :: send_req, recv_req
671 integer,
allocatable :: recv_buffer(:)
672 integer :: i, j, k, el_glb_idx, n_sides, n_nodes, facet, element, l
673 integer :: max_recv, ierr, src, dst, n_recv, recv_side, neigh_el
676 if (this%gdim .eq. 2)
then
689 el_glb_idx = i + this%offset_el
692 if (this%facet_neigh(j, i) .eq. 0)
then
693 if (n_nodes .eq. 2)
then
694 call this%elements(i)%e%facet_id(edge, j)
695 call buffer%push(el_glb_idx)
696 call buffer%push(facet)
698 call buffer%push(edge%x(k))
701 call this%elements(i)%e%facet_id(face, j)
702 call buffer%push(el_glb_idx)
703 call buffer%push(facet)
705 call buffer%push(face%x(k))
713 call mpi_allreduce(buffer%size(), max_recv, 1, &
716 allocate(recv_buffer(max_recv))
718 do i = 1,
size(this%neigh_order)
722 if (this%neigh(src))
then
723 call mpi_irecv(recv_buffer, max_recv, mpi_integer, &
727 if (this%neigh(dst))
then
728 call mpi_isend(buffer%array(), buffer%size(), mpi_integer, &
732 if (this%neigh(src))
then
733 call mpi_wait(recv_req, status, ierr)
734 call mpi_get_count(status, mpi_integer, n_recv, ierr)
736 select type (fmp => this%facet_map)
738 do j = 1, n_recv, n_nodes + 2
739 neigh_el = recv_buffer(j)
740 recv_side = recv_buffer(j+1)
742 edge = (/ recv_buffer(j+2), recv_buffer(j+3) /)
744 facet_data = (/ 0, 0 /)
746 if (fmp%get(edge, facet_data) .eq. 0)
then
747 element = facet_data%x(1) - this%offset_el
750 call this%elements(
element)%e%facet_id(edge2, l)
751 if(edge2 .eq. edge)
then
756 this%facet_neigh(facet,
element) = -neigh_el
757 facet_data%x(2) = -neigh_el
760 call fmp%set(edge, facet_data)
764 if (this%hte%get(edge, facet) .eq. 0)
then
774 do j = 1, n_recv, n_nodes + 2
775 neigh_el = recv_buffer(j)
776 recv_side = recv_buffer(j+1)
778 face%x = (/ recv_buffer(j+2), recv_buffer(j+3), &
779 recv_buffer(j+4), recv_buffer(j+5) /)
782 facet_data%x = (/ 0, 0 /)
785 if (fmp%get(face, facet_data) .eq. 0)
then
787 element = facet_data%x(1) - this%offset_el
789 call this%elements(
element)%e%facet_id(face2, l)
790 if(face2 .eq. face)
then
795 this%facet_neigh(facet,
element) = -neigh_el
796 facet_data%x(2) = -neigh_el
799 call fmp%set(face, facet_data)
803 if (this%htf%get(face, facet) .eq. 0)
then
816 if (this%neigh(dst))
then
817 call mpi_wait(send_req, mpi_status_ignore, ierr)
823 deallocate(recv_buffer)
831 type(
mesh_t),
intent(inout) :: this
833 type(mpi_status) :: status
834 integer,
allocatable :: recv_buffer(:)
836 integer :: max_recv, ierr, src, dst, n_recv, neigh_el
837 integer :: pt_glb_idx, pt_loc_idx, num_neigh
838 integer,
contiguous,
pointer :: neighs(:)
841 call send_buffer%init(this%mpts * 2)
846 pt_glb_idx = this%points(i)%id()
847 num_neigh = this%point_neigh(i)%size()
848 call send_buffer%push(pt_glb_idx)
849 call send_buffer%push(num_neigh)
851 neighs => this%point_neigh(i)%array()
852 do j = 1, this%point_neigh(i)%size()
853 call send_buffer%push(neighs(j))
857 call mpi_allreduce(send_buffer%size(), max_recv, 1, &
859 allocate(recv_buffer(max_recv))
865 call mpi_sendrecv(send_buffer%array(), send_buffer%size(), &
866 mpi_integer, dst, 0, recv_buffer, max_recv, mpi_integer, src, 0, &
869 call mpi_get_count(status, mpi_integer, n_recv, ierr)
872 do while (j .le. n_recv)
873 pt_glb_idx = recv_buffer(j)
874 num_neigh = recv_buffer(j + 1)
876 pt_loc_idx = this%have_point_glb_idx(pt_glb_idx)
877 if (pt_loc_idx .gt. 0)
then
878 this%neigh(src) = .true.
880 neigh_el = -recv_buffer(j + 1 + k)
881 call this%point_neigh(pt_loc_idx)%push(neigh_el)
885 j = j + (2 + num_neigh)
890 deallocate(recv_buffer)
891 call send_buffer%free()
899 type(
mesh_t),
target,
intent(inout) :: this
902 type(
uset_i8_t),
target :: edge_idx, ghost, owner
905 type(mpi_status) :: status
906 type(mpi_request) :: send_req, recv_req
907 integer,
contiguous,
pointer :: p1(:), p2(:), ns_id(:)
908 integer :: i, j, id, ierr, num_edge_glb, edge_offset, num_edge_loc
909 integer :: k, l , shared_offset, glb_nshared, n_glb_id
910 integer(kind=i8) :: C, glb_max, glb_id
911 integer(kind=i8),
pointer :: glb_ptr
912 integer(kind=i8),
allocatable :: recv_buff(:)
913 logical :: shared_edge
915 integer :: max_recv, src, dst, n_recv
919 allocate(this%ddata%local_to_global_edge(this%meds))
921 call edge_idx%init(this%hte%num_entries())
922 call send_buff%init(this%hte%num_entries())
923 call owner%init(this%hte%num_entries())
925 call glb_to_loc%init(32, i)
933 num_edge_glb = 2* this%meds
934 call mpi_allreduce(mpi_in_place, num_edge_glb, 1, &
937 glb_max = int(num_edge_glb,
i8)
939 call non_shared_edges%init(this%hte%num_entries())
941 call it%init(this%hte)
946 k = this%have_point_glb_idx(edge%x(1))
947 l = this%have_point_glb_idx(edge%x(2))
948 p1 => this%point_neigh(k)%array()
949 p2 => this%point_neigh(l)%array()
951 shared_edge = .false.
954 do i = 1, this%point_neigh(k)%size()
955 do j = 1, this%point_neigh(l)%size()
956 if ((p1(i) .eq. p2(j)) .and. &
957 (p1(i) .lt. 0) .and. (p2(j) .lt. 0))
then
967 if (shared_edge)
then
968 glb_id = ((int(edge%x(1),
i8)) + int(edge%x(2),
i8)*c) + glb_max
969 call glb_to_loc%set(glb_id, id)
970 call edge_idx%add(glb_id)
971 call owner%add(glb_id)
972 call send_buff%push(glb_id)
974 call non_shared_edges%push(id)
980 num_edge_loc = non_shared_edges%size()
981 call mpi_exscan(num_edge_loc, edge_offset, 1, &
983 edge_offset = edge_offset + 1
986 ns_id => non_shared_edges%array()
987 do i = 1, non_shared_edges%size()
989 edge_offset = edge_offset + 1
997 call mpi_allreduce(send_buff%size(), max_recv, 1, &
1000 call ghost%init(send_buff%size())
1002 allocate(recv_buff(max_recv))
1004 do i = 1,
size(this%neigh_order)
1008 if (this%neigh(src))
then
1009 call mpi_irecv(recv_buff, max_recv, mpi_integer8, &
1013 if (this%neigh(dst))
then
1017 select type(sbarray=>send_buff%data)
1018 type is (
integer(i8))
1019 call mpi_isend(sbarray, send_buff%size(), mpi_integer8, &
1024 if (this%neigh(src))
then
1025 call mpi_wait(recv_req, status, ierr)
1026 call mpi_get_count(status, mpi_integer8, n_recv, ierr)
1029 if ((edge_idx%element(recv_buff(j))) .and. (src .lt.
pe_rank))
then
1030 call ghost%add(recv_buff(j))
1031 call owner%remove(recv_buff(j))
1036 if (this%neigh(dst))
then
1037 call mpi_wait(send_req, mpi_status_ignore, ierr)
1043 glb_nshared = num_edge_loc
1044 call mpi_allreduce(mpi_in_place, glb_nshared, 1, &
1048 call mpi_exscan(owner%size(), shared_offset, 1, &
1050 shared_offset = shared_offset + glb_nshared + 1
1053 call send_buff%clear()
1054 call owner%iter_init()
1055 do while (owner%iter_next())
1056 glb_ptr => owner%iter_value()
1057 if (glb_to_loc%get(glb_ptr, id) .eq. 0)
then
1061 call send_buff%push(glb_ptr)
1062 glb_id = int(shared_offset, i8)
1063 call send_buff%push(glb_id)
1065 shared_offset = shared_offset + 1
1074 this%glb_meds = shared_offset -1
1075 call mpi_allreduce(mpi_in_place, this%glb_meds, 1, &
1082 call mpi_allreduce(send_buff%size(), max_recv, 1, &
1085 deallocate(recv_buff)
1086 allocate(recv_buff(max_recv))
1089 do i = 1,
size(this%neigh_order)
1093 if (this%neigh(src))
then
1094 call mpi_irecv(recv_buff, max_recv, mpi_integer8, &
1098 if (this%neigh(dst))
then
1102 select type(sbarray=>send_buff%data)
1103 type is (
integer(i8))
1104 call mpi_isend(sbarray, send_buff%size(), mpi_integer8, &
1109 if (this%neigh(src))
then
1110 call mpi_wait(recv_req, status, ierr)
1111 call mpi_get_count(status, mpi_integer8, n_recv, ierr)
1114 if (ghost%element(recv_buff(j)))
then
1115 if (glb_to_loc%get(recv_buff(j), id) .eq. 0)
then
1116 n_glb_id = int(recv_buff(j + 1 ), 4)
1126 if (this%neigh(dst))
then
1127 call mpi_wait(send_req, mpi_status_ignore, ierr)
1131 deallocate(recv_buff)
1132 call glb_to_loc%free()
1133 call send_buff%free()
1134 call edge_idx%free()
1135 call non_shared_edges%free()
1143 type(
mesh_t),
target,
intent(inout) :: this
1156 type(mpi_status) :: status
1157 type(mpi_request) :: send_req, recv_req
1158 integer,
allocatable :: recv_buff(:)
1159 integer :: non_shared_facets, shared_facets, facet_offset
1160 integer :: id, glb_nshared, shared_offset, owned_facets
1161 integer :: i, j, ierr, max_recv, src, dst, n_recv
1163 shared_facets = this%ddata%shared_facet%size()
1166 if (this%gdim .eq. 2)
then
1167 allocate(this%ddata%local_to_global_facet(this%meds))
1168 call edge_owner%init(this%meds)
1169 call edge_ghost%init(64, i)
1170 non_shared_facets = this%hte%num_entries() - shared_facets
1172 allocate(this%ddata%local_to_global_facet(this%mfcs))
1173 call face_owner%init(this%mfcs)
1174 call face_ghost%init(64, i)
1175 non_shared_facets = this%htf%num_entries() - shared_facets
1181 call mpi_exscan(non_shared_facets, facet_offset, 1, &
1183 facet_offset = facet_offset + 1
1186 if (this%gdim .eq. 2)
then
1187 call edge_it%init(this%hte)
1188 do while (edge_it%next())
1189 call edge_it%data(id)
1190 edge => edge_it%key()
1191 if (.not. this%ddata%shared_facet%element(id))
then
1194 facet_offset = facet_offset + 1
1196 select type(fmp => this%facet_map)
1198 if (fmp%get(edge, facet_data) .eq. 0)
then
1199 if (facet_data%x(2) .lt. 0)
then
1200 if (abs(facet_data%x(2)) .lt. (this%offset_el + 1))
then
1201 call edge_ghost%set(edge, id)
1203 call edge_owner%push(edge)
1212 owned_facets = edge_owner%size()
1214 call face_it%init(this%htf)
1215 do while (face_it%next())
1216 call face_it%data(id)
1217 face => face_it%key()
1218 if (.not. this%ddata%shared_facet%element(id))
then
1221 facet_offset = facet_offset + 1
1223 select type(fmp => this%facet_map)
1225 if (fmp%get(face, facet_data) .eq. 0)
then
1226 if (facet_data%x(2) .lt. 0)
then
1227 if (abs(facet_data%x(2)) .lt. (this%offset_el + 1))
then
1228 call face_ghost%set(face, id)
1230 call face_owner%push(face)
1239 owned_facets = face_owner%size()
1243 glb_nshared = non_shared_facets
1244 call mpi_allreduce(mpi_in_place, glb_nshared, 1, &
1248 call mpi_exscan(owned_facets, shared_offset, 1, &
1250 shared_offset = shared_offset + glb_nshared + 1
1252 if (this%gdim .eq. 2)
then
1254 if (owned_facets .gt. 32)
then
1255 call send_buff%init(owned_facets)
1257 call send_buff%init()
1260 ed => edge_owner%array()
1261 do i = 1, edge_owner%size()
1262 if (this%hte%get(ed(i), id) .eq. 0)
then
1269 call send_buff%push(ed(i)%x(j))
1271 call send_buff%push(shared_offset)
1273 shared_offset = shared_offset + 1
1279 if (owned_facets .gt. 32)
then
1280 call send_buff%init(owned_facets)
1282 call send_buff%init()
1285 fd => face_owner%array()
1286 do i = 1, face_owner%size()
1287 if (this%htf%get(fd(i), id) .eq. 0)
then
1294 call send_buff%push(fd(i)%x(j))
1296 call send_buff%push(shared_offset)
1298 shared_offset = shared_offset + 1
1307 this%glb_mfcs = shared_offset - 1
1308 call mpi_allreduce(mpi_in_place, this%glb_mfcs, 1, &
1315 call mpi_allreduce(send_buff%size(), max_recv, 1, &
1318 allocate(recv_buff(max_recv))
1321 do i = 1,
size(this%neigh_order)
1325 if (this%neigh(src))
then
1326 call mpi_irecv(recv_buff, max_recv, mpi_integer, &
1330 if (this%neigh(dst))
then
1331 call mpi_isend(send_buff%array(), send_buff%size(), mpi_integer, &
1335 if (this%neigh(src))
then
1336 call mpi_wait(recv_req, status, ierr)
1337 call mpi_get_count(status, mpi_integer, n_recv, ierr)
1339 if (this%gdim .eq. 2)
then
1342 recv_edge = (/recv_buff(j), recv_buff(j+1)/)
1345 if (edge_ghost%get(recv_edge, id) .eq. 0)
then
1353 recv_face = (/recv_buff(j), recv_buff(j+1), &
1354 recv_buff(j+2), recv_buff(j+3) /)
1357 if (face_ghost%get(recv_face, id) .eq. 0)
then
1365 if (this%neigh(dst))
then
1366 call mpi_wait(send_req, mpi_status_ignore, ierr)
1371 if (this%gdim .eq. 2)
then
1372 call edge_owner%free()
1373 call edge_ghost%free()
1375 call face_owner%free()
1376 call face_ghost%free()
1379 call send_buff%free()
1380 deallocate(recv_buff)
1387 class(
mesh_t),
target,
intent(inout) :: this
1388 integer,
value :: el
1389 type(
point_t),
target,
intent(inout) :: p1, p2, p3, p4
1391 integer :: p(4), el_glb_idx, i, p_local_idx
1395 this%lconn = .false.
1398 this%lnumr = .false.
1400 call this%add_point(p1, p(1))
1401 call this%add_point(p2, p(2))
1402 call this%add_point(p3, p(3))
1403 call this%add_point(p4, p(4))
1405 ep => this%elements(el)%e
1406 el_glb_idx = el + this%offset_el
1410 call ep%init(el_glb_idx, &
1411 this%points(p(1)), this%points(p(2)), &
1412 this%points(p(3)), this%points(p(4)))
1423 class(
mesh_t),
target,
intent(inout) :: this
1424 integer,
value :: el
1425 type(
point_t),
target,
intent(inout) :: p1, p2, p3, p4, p5, p6, p7, p8
1427 integer :: p(8), el_glb_idx, i, p_local_idx
1432 this%lconn = .false.
1435 this%lnumr = .false.
1437 call this%add_point(p1, p(1))
1438 call this%add_point(p2, p(2))
1439 call this%add_point(p3, p(3))
1440 call this%add_point(p4, p(4))
1441 call this%add_point(p5, p(5))
1442 call this%add_point(p6, p(6))
1443 call this%add_point(p7, p(7))
1444 call this%add_point(p8, p(8))
1446 ep => this%elements(el)%e
1447 el_glb_idx = el + this%offset_el
1450 call ep%init(el_glb_idx, &
1451 this%points(p(1)), this%points(p(2)), &
1452 this%points(p(3)), this%points(p(4)), &
1453 this%points(p(5)), this%points(p(6)), &
1454 this%points(p(7)), this%points(p(8)))
1463 class(
mesh_t),
intent(inout) :: this
1464 type(
point_t),
intent(inout) :: p
1465 integer,
intent(inout) :: idx
1470 this%max_pts_id =
max(this%max_pts_id, tmp)
1472 if (tmp .le. 0)
then
1476 if (this%htp%get(tmp, idx) .gt. 0)
then
1477 this%mpts = this%mpts + 1
1478 call this%htp%set(tmp, this%mpts)
1479 this%points(this%mpts) = p
1487 class(
mesh_t),
intent(inout) :: this
1491 if (this%htf%get(f, idx) .gt. 0)
then
1492 this%mfcs = this%mfcs + 1
1493 call this%htf%set(f, this%mfcs)
1500 class(
mesh_t),
intent(inout) :: this
1504 if (this%hte%get(e, idx) .gt. 0)
then
1505 this%meds = this%meds + 1
1506 call this%hte%set(e, this%meds)
1513 class(
mesh_t),
intent(inout) :: this
1514 integer,
intent(inout) :: f
1515 integer,
intent(inout) :: e
1517 if (e .gt. this%nelv)
then
1521 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1522 (this%gdim .eq. 3 .and. f .gt. 6))
then
1525 this%facet_type(f, e) = 2
1526 call this%wall%add_facet(f, e)
1532 class(
mesh_t),
intent(inout) :: this
1533 integer,
intent(in) :: e
1534 real(kind=
dp),
dimension(5,12),
intent(in) :: curve_data
1535 integer,
dimension(12),
intent(in) :: curve_type
1537 if (e .gt. this%nelv)
then
1540 if ((this%gdim .eq. 2 .and. sum(curve_type(5:8)) .gt. 0) )
then
1543 call this%curve%add_element(e, curve_data, curve_type)
1550 class(
mesh_t),
intent(inout) :: this
1551 integer,
intent(in) :: f
1552 integer,
intent(in) :: e
1554 if (e .gt. this%nelv)
then
1558 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1559 (this%gdim .eq. 3 .and. f .gt. 6))
then
1562 this%facet_type(f, e) = 2
1563 call this%inlet%add_facet(f, e)
1569 class(
mesh_t),
intent(inout) :: this
1570 integer,
intent(in) :: f
1571 integer,
intent(in) :: e
1572 integer,
intent(in) :: label
1574 if (e .gt. this%nelv)
then
1578 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1579 (this%gdim .eq. 3 .and. f .gt. 6))
then
1582 call this%labeled_zones(label)%add_facet(f, e)
1583 this%facet_type(f,e) = -label
1590 class(
mesh_t),
intent(inout) :: this
1591 integer,
intent(inout) :: f
1592 integer,
intent(inout) :: e
1594 if (e .gt. this%nelv)
then
1598 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1599 (this%gdim .eq. 3 .and. f .gt. 6))
then
1602 this%facet_type(f, e) = 1
1603 call this%outlet_normal%add_facet(f, e)
1610 class(
mesh_t),
intent(inout) :: this
1611 integer,
intent(inout) :: f
1612 integer,
intent(inout) :: e
1614 if (e .gt. this%nelv)
then
1618 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1619 (this%gdim .eq. 3 .and. f .gt. 6))
then
1622 this%facet_type(f, e) = 1
1623 call this%outlet%add_facet(f, e)
1629 class(
mesh_t),
intent(inout) :: this
1630 integer,
intent(inout) :: f
1631 integer,
intent(inout) :: e
1633 if (e .gt. this%nelv)
then
1637 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1638 (this%gdim .eq. 3 .and. f .gt. 6))
then
1641 this%facet_type(f, e) = 2
1642 call this%sympln%add_facet(f, e)
1648 class(
mesh_t),
intent(inout) :: this
1649 integer,
intent(in) :: f
1650 integer,
intent(in) :: e
1651 integer,
intent(in) :: pf
1652 integer,
intent(in) :: pe
1653 integer,
intent(inout) :: pids(4)
1654 integer,
dimension(4) :: org_ids
1656 call this%get_facet_ids(f, e, org_ids)
1657 call this%periodic%add_periodic_facet(f, e, pf, pe, pids, org_ids)
1662 class(
mesh_t),
intent(inout) :: this
1663 integer,
intent(in) :: f
1664 integer,
intent(in) :: e
1665 integer,
intent(inout) :: pids(4)
1671 select type(ele => this%elements(e)%e)
1673 call ele%facet_order(t,f)
1676 call ele%facet_order(t2,f)
1686 class(
mesh_t),
intent(inout) :: this
1692 integer :: org_ids(4), pids(4)
1694 integer,
dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1701 integer,
dimension(2, 4) :: edge_nodes = reshape((/1,3,&
1707 do i = 1, this%periodic%size
1708 e = this%periodic%facet_el(i)%x(2)
1709 f = this%periodic%facet_el(i)%x(1)
1710 pe = this%periodic%p_facet_el(i)%x(2)
1711 pf = this%periodic%p_facet_el(i)%x(1)
1712 pids = this%periodic%p_ids(i)%x
1713 call this%get_facet_ids(f, e, pids)
1714 this%periodic%p_ids(i)%x = pids
1716 do i = 1, this%periodic%size
1717 e = this%periodic%facet_el(i)%x(2)
1718 f = this%periodic%facet_el(i)%x(1)
1719 org_ids = this%periodic%org_ids(i)%x
1720 select type(ele => this%elements(e)%e)
1723 pi => ele%pts(face_nodes(j,f))%p
1724 call pi%set_id(org_ids(j))
1728 pi => ele%pts(edge_nodes(j,f))%p
1729 call pi%set_id(org_ids(j))
1737 class(
mesh_t),
intent(inout) :: this
1738 integer,
intent(in) :: f
1739 integer,
intent(in) :: e
1740 integer,
intent(in) :: pf
1741 integer,
intent(in) :: pe
1742 type(
point_t),
pointer :: pi, pj
1743 real(kind=
dp) :: l(3)
1744 integer :: i, j, id, p_local_idx, match
1747 integer,
dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1754 integer,
dimension(2, 4) :: edge_nodes = reshape((/1,3,&
1760 select type(ele => this%elements(e)%e)
1762 select type(elp => this%elements(pe)%e)
1766 l = l + ele%pts(face_nodes(i,f))%p%x(1:3) - &
1767 elp%pts(face_nodes(i,pf))%p%x(1:3)
1771 pi => ele%pts(face_nodes(i,f))%p
1774 pj => elp%pts(face_nodes(j,pf))%p
1775 if (norm2(pi%x(1:3) - pj%x(1:3) - l) .lt. 1d-7)
then
1776 id = min(pi%id(), pj%id())
1779 p_local_idx = this%get_local(this%points(id))
1783 if ( match .gt. 1)
then
1784 call neko_error(
'Multiple matches when creating periodic ids')
1785 else if (match .eq. 0)
then
1786 call neko_error(
'Cannot find matching periodic point')
1791 select type(elp => this%elements(pe)%e)
1795 l = l + ele%pts(edge_nodes(i,f))%p%x(1:3) - &
1796 elp%pts(edge_nodes(i,pf))%p%x(1:3)
1800 pi => ele%pts(edge_nodes(i,f))%p
1802 pj => elp%pts(edge_nodes(j,pf))%p
1804 if (norm2(pi%x(1:3) - pj%x(1:3) - l) .lt. 1d-7)
then
1805 id = min(pi%id(), pj%id())
1808 p_local_idx = this%get_local(this%points(id))
1819 class(
mesh_t),
intent(inout) :: this
1820 integer,
intent(in) :: f
1821 integer,
intent(in) :: e
1822 integer,
intent(in) :: pf
1823 integer,
intent(in) :: pe
1824 integer,
intent(inout) :: pids(4)
1826 integer :: i, id, p_local_idx
1829 integer,
dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1836 select type(ele => this%elements(e)%e)
1839 pi => ele%pts(face_nodes(i,f))%p
1840 call pi%set_id(pids(i))
1841 call this%add_point(pi, id)
1842 p_local_idx = this%get_local(this%points(id))
1850 class(
mesh_t),
intent(inout) :: this
1851 type(
point_t),
intent(inout) :: p
1858 if (this%htp%get(tmp, local_id) .gt. 0)
then
1859 call neko_error(
'Invalid global id (local point)')
1867 class(
mesh_t),
intent(inout) :: this
1871 if (this%hte%get(e, local_id) .gt. 0)
then
1872 call neko_error(
'Invalid global id (local edge)')
1879 class(
mesh_t),
intent(inout) :: this
1883 if (this%htf%get(f, local_id) .gt. 0)
then
1884 call neko_error(
'Invalid global id (local facet)')
1891 class(
mesh_t),
intent(inout) :: this
1893 integer :: global_id
1895 global_id = this%get_local(e)
1897 if (this%gdim .eq. 2)
then
1899 global_id = this%ddata%local_to_global_facet(global_id)
1903 global_id = this%ddata%local_to_global_edge(global_id)
1911 class(
mesh_t),
intent(inout) :: this
1913 integer :: global_id
1915 global_id = this%get_local_facet(f)
1918 global_id = this%ddata%local_to_global_facet(global_id)
1928 class(
mesh_t),
intent(inout) :: this
1929 integer,
intent(inout) :: index
1932 if (this%htp%get(index, local_id) .eq. 1)
then
1941 class(
mesh_t),
intent(inout) :: this
1942 type(
point_t),
intent(inout) :: p
1943 integer :: local_index
1946 local_index = this%get_local(p)
1947 shared = this%ddata%shared_point%element(local_index)
1955 class(
mesh_t),
intent(inout) :: this
1957 integer :: local_index
1959 local_index = this%get_local(e)
1960 if (this%gdim .eq. 2)
then
1961 shared = this%ddata%shared_facet%element(local_index)
1963 shared = this%ddata%shared_edge%element(local_index)
1969 class(
mesh_t),
intent(inout) :: this
1971 integer :: local_index
1974 local_index = this%get_local(f)
1975 shared = this%ddata%shared_facet%element(local_index)
type(mpi_comm) neko_comm
MPI communicator.
integer pe_size
MPI size of communicator.
Defines a domain as a subset of facets in a mesh.
Defines practical data distributions.
subroutine, public distdata_set_shared_point(ddata, point)
Mark a point as shared.
subroutine, public distdata_set_shared_el_facet(ddata, element, side)
Mark an element's facet as shared.
subroutine, public distdata_set_local_to_global_facet(ddata, local, global)
Set local to global mapping (facets)
subroutine, public distdata_free(ddata)
Free a distdata type.
subroutine, public distdata_set_local_to_global_edge(ddata, local, global)
Set local to global mapping (edges)
subroutine, public distdata_init(ddata)
Initialise a distdata type.
subroutine, public distdata_set_shared_facet(ddata, facet)
Mark a facet as shared.
subroutine, public distdata_set_shared_edge(ddata, edge)
Mark an element's edge as shared.
Defines a zone as a subset of facets in a mesh.
Defines a hexahedron element.
integer, parameter, public neko_hex_npts
Number of points.
integer, parameter, public neko_hex_nfcs
Number of faces.
integer, parameter, public neko_hex_neds
Number of edges.
Implements a hash table ADT.
subroutine mesh_generate_flags(this)
subroutine mesh_generate_facet_numbering(this)
Generate a unique facet numbering.
subroutine mesh_init_dist(this, gdim, dist)
Initialise a mesh this based on a distribution dist.
logical function mesh_is_shared_point(this, p)
Check if a point is shared.
subroutine mesh_mark_wall_facet(this, f, e)
Mark facet f in element e as a wall.
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
subroutine mesh_mark_sympln_facet(this, f, e)
Mark facet f in element e as a symmetry plane.
subroutine mesh_add_face(this, f)
Add a unique face represented as a 4-tuple to the mesh.
integer function mesh_get_global_edge(this, e)
Return the global id of an edge e.
subroutine mesh_generate_edge_conn(this)
Generate element-element connectivity via edges both between internal and between PEs.
subroutine mesh_add_hex(this, el, p1, p2, p3, p4, p5, p6, p7, p8)
Add a hexahedral element to the mesh this.
subroutine mesh_add_edge(this, e)
Add a unique edge represented as a 2-tuple to the mesh.
subroutine mesh_add_point(this, p, idx)
Add a unique point to the mesh.
subroutine mesh_free(this)
Deallocate a mesh this.
subroutine mesh_mark_labeled_facet(this, f, e, label)
Mark facet f in element e with label.
integer, parameter, public neko_msh_max_zlbl_len
Max length of a zone label.
subroutine mesh_mark_periodic_facet(this, f, e, pf, pe, pids)
Mark facet f in element e as periodic with (pf, pe)
integer function mesh_get_local_edge(this, e)
Return the local id of an edge e.
subroutine mesh_mark_outlet_normal_facet(this, f, e)
Mark facet f in element e as an outlet normal.
subroutine mesh_all_deformed(this)
Set all elements as if they are deformed.
logical function mesh_is_shared_edge(this, e)
Check if an edge is shared.
subroutine mesh_mark_inlet_facet(this, f, e)
Mark facet f in element e as an inlet.
subroutine mesh_init_common(this)
subroutine mesh_get_facet_ids(this, f, e, pids)
Get original ids of periodic points.
subroutine mesh_create_periodic_ids(this, f, e, pf, pe)
Creates common ids for matching periodic points.
subroutine mesh_reset_periodic_ids(this)
Reset ids of periodic points to their original ids.
subroutine mesh_init_nelv(this, gdim, nelv)
Initialise a mesh this with nelv elements.
logical function mesh_is_shared_facet(this, f)
Check if a facet is shared.
integer function mesh_get_global_facet(this, f)
Return the local id of a face f.
integer function mesh_have_point_glb_idx(this, index)
Check if the mesh has a point given its global index.
subroutine mesh_generate_external_point_conn(this)
Generate element-element connectivity via points between PEs.
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.
subroutine mesh_finalize(this)
subroutine mesh_generate_external_facet_conn(this)
Generate element-element connectivity via facets between PEs.
subroutine mesh_mark_curve_element(this, e, curve_data, curve_type)
Mark element e as a curve element.
subroutine mesh_generate_conn(this)
Generate element-to-element connectivity.
subroutine mesh_mark_outlet_facet(this, f, e)
Mark facet f in element e as an outlet.
integer function mesh_get_local_facet(this, f)
Return the local id of a face f.
integer function mesh_get_local_point(this, p)
Return the local id of a point p.
subroutine mesh_add_quad(this, el, p1, p2, p3, p4)
Add a quadrilateral element to the mesh this.
integer, parameter, public i8
integer, parameter, public dp
integer, parameter, public rp
Global precision used in computations.
Defines a quadrilateral element.
integer, parameter, public neko_quad_neds
Number of edges.
integer, parameter, public neko_quad_npts
Number of points.
Implements a dynamic stack ADT.
Implements an unordered set ADT.
Defines wall boundary conditions.
Load-balanced linear distribution .
Base type for an element.
Integer based hash table.
Integer 2-tuple based hash table.
Integer 4-tuple based hash table.
Integer*8 based hash table.
Iterator for an integer based 2-tuple hash table.
Iterator for an integer based 4-tuple hash table.
Base type for a hash table.
A point in with coordinates .
Integer 2-tuple based stack.
Integer 4-tuple based stack.
Integer*8 based unordered set.