82 logical,
allocatable :: dfrmd_el(:)
88 integer,
allocatable :: facet_neigh(:,:)
97 logical,
allocatable :: neigh(:)
98 integer,
allocatable :: neigh_order(:)
100 integer(2),
allocatable :: facet_type(:,:)
111 logical :: lconn = .false.
112 logical :: ldist = .false.
113 logical :: lnumr = .false.
114 logical :: lgenc = .true.
118 procedure(
mesh_deform), pass(msh),
pointer :: apply_deform => null()
142 procedure, pass(this) :: mark_outlet_normal_facet => &
154 generic :: init => init_nelv, init_dist
156 generic :: add_element => add_quad, add_hex
159 generic :: get_local => get_local_point, get_local_edge, get_local_facet
162 generic :: get_global => get_global_edge, get_global_facet
164 generic :: is_shared => is_shared_point, is_shared_edge, is_shared_facet
172 integer,
intent(in) :: lx, ly, lz
173 real(kind=
rp),
intent(inout) :: x(lx, ly, lz, msh%nelv)
174 real(kind=
rp),
intent(inout) :: y(lx, ly, lz, msh%nelv)
175 real(kind=
rp),
intent(inout) :: z(lx, ly, lz, msh%nelv)
185 class(
mesh_t),
intent(inout) :: this
186 integer,
intent(in) :: gdim
187 integer,
intent(in) :: nelv
189 character(len=LOG_SIZE) :: log_buf
196 if (this%nelv < 1)
then
197 write(log_buf,
'(A,I0,A)')
'MPI rank ',
pe_rank,
' has zero elements'
201 call mpi_allreduce(this%nelv, this%glb_nelv, 1, &
205 call mpi_exscan(this%nelv, this%offset_el, 1, &
214 class(
mesh_t),
intent(inout) :: this
215 integer,
intent(in) :: gdim
217 character(len=LOG_SIZE) :: log_buf
221 this%nelv = dist%num_local()
222 if (this%nelv < 1)
then
223 write(log_buf,
'(A,I0,A)')
'MPI rank ',
pe_rank,
' has zero elements'
226 this%glb_nelv = dist%num_global()
227 this%offset_el = dist%start_idx()
235 type(
mesh_t),
intent(inout) :: this
241 allocate(this%elements(this%nelv))
242 allocate(this%dfrmd_el(this%nelv))
243 if (this%gdim .eq. 3)
then
245 allocate(
hex_t::this%elements(i)%e)
251 select type (fmp => this%facet_map)
253 call fmp%init(this%nelv, facet_data)
261 else if (this%gdim .eq. 2)
then
263 allocate(
quad_t::this%elements(i)%e)
268 select type (fmp => this%facet_map)
270 call fmp%init(this%nelv, facet_data)
282 allocate(this%points(this%npts*this%nelv))
287 allocate(this%point_neigh(this%gdim*this%npts*this%nelv))
288 do i = 1, this%gdim*this%npts*this%nelv
289 call this%point_neigh(i)%init()
293 allocate(this%facet_type(2 * this%gdim, this%nelv))
296 call this%htp%init(this%npts*this%nelv, i)
298 call this%wall%init(this%nelv)
299 call this%inlet%init(this%nelv)
300 call this%outlet%init(this%nelv)
301 call this%outlet_normal%init(this%nelv)
302 call this%sympln%init(this%nelv)
303 call this%periodic%init(this%nelv)
307 call this%labeled_zones(i)%init(this%nelv)
310 call this%curve%init(this%nelv)
314 allocate(this%neigh(0:
pe_size-1))
325 class(
mesh_t),
intent(inout) :: this
334 if (
allocated(this%dfrmd_el))
then
335 deallocate(this%dfrmd_el)
338 if (
allocated(this%elements))
then
340 call this%elements(i)%e%free()
341 deallocate(this%elements(i)%e)
343 deallocate(this%elements)
346 if (
allocated(this%facet_map))
then
347 select type (fmp => this%facet_map)
353 deallocate(this%facet_map)
356 if (
allocated(this%facet_neigh))
then
357 deallocate(this%facet_neigh)
360 if (
allocated(this%point_neigh))
then
361 do i = 1, this%gdim * this%npts * this%nelv
362 call this%point_neigh(i)%free()
364 deallocate(this%point_neigh)
367 if (
allocated(this%facet_type))
then
368 deallocate(this%facet_type)
370 if (
allocated(this%labeled_zones))
then
372 call this%labeled_zones(i)%free()
374 deallocate(this%labeled_zones)
377 if (
allocated(this%neigh))
then
378 deallocate(this%neigh)
381 if (
allocated(this%neigh_order))
then
382 deallocate(this%neigh_order)
385 if (
allocated(this%points))
then
386 deallocate(this%points)
389 call this%wall%free()
390 call this%inlet%free()
391 call this%outlet%free()
392 call this%outlet_normal%free()
393 call this%sympln%free()
394 call this%periodic%free()
403 class(
mesh_t),
target,
intent(inout) :: this
409 call this%wall%finalize()
410 call this%inlet%finalize()
411 call this%outlet%finalize()
412 call this%outlet_normal%finalize()
413 call this%sympln%finalize()
414 call this%periodic%finalize()
416 call this%labeled_zones(i)%finalize()
418 call this%curve%finalize()
423 type(
mesh_t),
intent(inout) :: this
424 real(kind=
dp) :: u(3), v(3), w(3), temp
428 if (this%gdim .eq. 2)
then
429 this%dfrmd_el(e) = .false.
430 u = this%elements(e)%e%pts(2)%p%x - this%elements(e)%e%pts(1)%p%x
431 v = this%elements(e)%e%pts(3)%p%x - this%elements(e)%e%pts(1)%p%x
432 temp = u(1)*v(1) + u(2)*v(2)
433 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
435 this%dfrmd_el(e) = .false.
436 u = this%elements(e)%e%pts(2)%p%x - this%elements(e)%e%pts(1)%p%x
437 v = this%elements(e)%e%pts(3)%p%x - this%elements(e)%e%pts(1)%p%x
438 w = this%elements(e)%e%pts(5)%p%x - this%elements(e)%e%pts(1)%p%x
439 temp = u(1)*v(1) + u(2)*v(2) + u(3)*v(3)
440 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
441 temp = u(1)*w(1) + u(2)*w(2) + u(3)*w(3)
442 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
443 u = this%elements(e)%e%pts(7)%p%x - this%elements(e)%e%pts(8)%p%x
444 v = this%elements(e)%e%pts(6)%p%x - this%elements(e)%e%pts(8)%p%x
445 w = this%elements(e)%e%pts(4)%p%x - this%elements(e)%e%pts(8)%p%x
446 temp = u(1)*v(1) + u(2)*v(2) + u(3)*v(3)
447 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
448 temp = u(1)*w(1) + u(2)*w(2) + u(3)*w(3)
449 if(.not.
abscmp(temp, 0d0)) this%dfrmd_el(e) = .true.
456 class(
mesh_t),
intent(inout) :: this
457 this%dfrmd_el = .true.
462 class(
mesh_t),
target,
intent(inout) :: this
470 integer :: p_local_idx
472 integer :: i, j, k, ierr, el_glb_idx, n_sides, n_nodes, src, dst
474 if (this%lconn)
return
476 if (.not. this%lgenc)
return
480 ep => this%elements(el)%e
485 call this%add_point(ep%pts(i)%p, id)
486 p_local_idx = this%get_local(this%points(id))
489 call this%point_neigh(p_local_idx)%push(id)
492 call ep%facet_id(f, i)
493 call this%add_face(f)
497 call ep%edge_id(e, i)
498 call this%add_edge(e)
503 call this%add_point(ep%pts(i)%p, id)
504 p_local_idx = this%get_local(this%points(id))
507 call this%point_neigh(p_local_idx)%push(id)
511 call ep%facet_id(e, i)
512 call this%add_edge(e)
518 if (this%gdim .eq. 2)
then
527 call mpi_allreduce(this%max_pts_id, this%glb_mpts, 1, &
537 select type (fmp => this%facet_map)
541 el_glb_idx = i + this%offset_el
543 call this%elements(i)%e%facet_id(edge, j)
546 facet_data%x = (/ 0, 0/)
549 if (fmp%get(edge, facet_data) .eq. 0)
then
551 if (facet_data%x(1) .eq. el_glb_idx )
then
552 this%facet_neigh(j, i) = facet_data%x(2)
553 else if( facet_data%x(2) .eq. el_glb_idx)
then
554 this%facet_neigh(j, i) = facet_data%x(1)
556 else if(facet_data%x(1) .gt. el_glb_idx)
then
557 facet_data%x(2) = facet_data%x(1)
558 facet_data%x(1) = el_glb_idx
559 this%facet_neigh(j, i) = facet_data%x(2)
560 call fmp%set(edge, facet_data)
561 else if(facet_data%x(1) .lt. el_glb_idx)
then
562 facet_data%x(2) = el_glb_idx
563 this%facet_neigh(j, i) = facet_data%x(1)
564 call fmp%set(edge, facet_data)
567 facet_data%x(1) = el_glb_idx
568 this%facet_neigh(j, i) = facet_data%x(2)
569 call fmp%set(edge, facet_data)
578 el_glb_idx = i + this%offset_el
580 call this%elements(i)%e%facet_id(face, j)
582 facet_data%x = (/ 0, 0/)
585 if (fmp%get(face, facet_data) .eq. 0)
then
587 if (facet_data%x(1) .eq. el_glb_idx )
then
588 this%facet_neigh(j, i) = facet_data%x(2)
589 call this%elements(i)%e%facet_id(face_comp, &
591 if (face_comp .eq. face)
then
592 facet_data%x(2) = el_glb_idx
593 this%facet_neigh(j, i) = facet_data%x(1)
594 call fmp%set(face, facet_data)
596 else if( facet_data%x(2) .eq. el_glb_idx)
then
597 this%facet_neigh(j, i) = facet_data%x(1)
599 else if(facet_data%x(1) .gt. el_glb_idx)
then
600 facet_data%x(2) = facet_data%x(1)
601 facet_data%x(1) = el_glb_idx
602 this%facet_neigh(j, i) = facet_data%x(2)
603 call fmp%set(face, facet_data)
604 else if(facet_data%x(1) .lt. el_glb_idx)
then
605 facet_data%x(2) = el_glb_idx
606 this%facet_neigh(j, i) = facet_data%x(1)
607 call fmp%set(face, facet_data)
610 facet_data%x(1) = el_glb_idx
611 this%facet_neigh(j, i) = 0
612 call fmp%set(face, facet_data)
637 if (this%neigh(src) .or. this%neigh(dst))
then
639 call neigh_order%push(j)
643 allocate(this%neigh_order(neigh_order%size()))
644 select type(order => neigh_order%data)
646 do i = 1, neigh_order%size()
647 this%neigh_order(i) = order(i)
650 call neigh_order%free()
654 allocate(this%neigh_order(1))
663 if (this%gdim .eq. 3)
then
676 type(
mesh_t),
intent(inout) :: this
681 type(mpi_status) :: status
682 type(mpi_request) :: send_req, recv_req
683 integer,
allocatable :: recv_buffer(:)
684 integer :: i, j, k, el_glb_idx, n_sides, n_nodes, facet, element, l
685 integer :: max_recv, ierr, src, dst, n_recv, recv_side, neigh_el
688 if (this%gdim .eq. 2)
then
701 el_glb_idx = i + this%offset_el
704 if (this%facet_neigh(j, i) .eq. 0)
then
705 if (n_nodes .eq. 2)
then
706 call this%elements(i)%e%facet_id(edge, j)
707 call buffer%push(el_glb_idx)
708 call buffer%push(facet)
710 call buffer%push(edge%x(k))
713 call this%elements(i)%e%facet_id(face, j)
714 call buffer%push(el_glb_idx)
715 call buffer%push(facet)
717 call buffer%push(face%x(k))
725 call mpi_allreduce(buffer%size(), max_recv, 1, &
728 allocate(recv_buffer(max_recv))
730 do i = 1,
size(this%neigh_order)
734 if (this%neigh(src))
then
735 call mpi_irecv(recv_buffer, max_recv, mpi_integer, &
739 if (this%neigh(dst))
then
740 call mpi_isend(buffer%array(), buffer%size(), mpi_integer, &
744 if (this%neigh(src))
then
745 call mpi_wait(recv_req, status, ierr)
746 call mpi_get_count(status, mpi_integer, n_recv, ierr)
748 select type (fmp => this%facet_map)
750 do j = 1, n_recv, n_nodes + 2
751 neigh_el = recv_buffer(j)
752 recv_side = recv_buffer(j+1)
754 edge = (/ recv_buffer(j+2), recv_buffer(j+3) /)
756 facet_data = (/ 0, 0 /)
758 if (fmp%get(edge, facet_data) .eq. 0)
then
759 element = facet_data%x(1) - this%offset_el
762 call this%elements(
element)%e%facet_id(edge2, l)
763 if(edge2 .eq. edge)
then
768 this%facet_neigh(facet,
element) = -neigh_el
769 facet_data%x(2) = -neigh_el
772 call fmp%set(edge, facet_data)
776 if (this%hte%get(edge, facet) .eq. 0)
then
786 do j = 1, n_recv, n_nodes + 2
787 neigh_el = recv_buffer(j)
788 recv_side = recv_buffer(j+1)
790 face%x = (/ recv_buffer(j+2), recv_buffer(j+3), &
791 recv_buffer(j+4), recv_buffer(j+5) /)
794 facet_data%x = (/ 0, 0 /)
797 if (fmp%get(face, facet_data) .eq. 0)
then
799 element = facet_data%x(1) - this%offset_el
801 call this%elements(
element)%e%facet_id(face2, l)
802 if(face2 .eq. face)
then
807 this%facet_neigh(facet,
element) = -neigh_el
808 facet_data%x(2) = -neigh_el
811 call fmp%set(face, facet_data)
815 if (this%htf%get(face, facet) .eq. 0)
then
828 if (this%neigh(dst))
then
829 call mpi_wait(send_req, mpi_status_ignore, ierr)
835 deallocate(recv_buffer)
843 type(
mesh_t),
intent(inout) :: this
845 type(mpi_status) :: status
846 integer,
allocatable :: recv_buffer(:)
848 integer :: max_recv, ierr, src, dst, n_recv, neigh_el
849 integer :: pt_glb_idx, pt_loc_idx, num_neigh
850 integer,
contiguous,
pointer :: neighs(:)
853 call send_buffer%init(this%mpts * 2)
858 pt_glb_idx = this%points(i)%id()
859 num_neigh = this%point_neigh(i)%size()
860 call send_buffer%push(pt_glb_idx)
861 call send_buffer%push(num_neigh)
863 neighs => this%point_neigh(i)%array()
864 do j = 1, this%point_neigh(i)%size()
865 call send_buffer%push(neighs(j))
869 call mpi_allreduce(send_buffer%size(), max_recv, 1, &
871 allocate(recv_buffer(max_recv))
877 call mpi_sendrecv(send_buffer%array(), send_buffer%size(), &
878 mpi_integer, dst, 0, recv_buffer, max_recv, mpi_integer, src, 0, &
881 call mpi_get_count(status, mpi_integer, n_recv, ierr)
884 do while (j .le. n_recv)
885 pt_glb_idx = recv_buffer(j)
886 num_neigh = recv_buffer(j + 1)
888 pt_loc_idx = this%have_point_glb_idx(pt_glb_idx)
889 if (pt_loc_idx .gt. 0)
then
890 this%neigh(src) = .true.
892 neigh_el = -recv_buffer(j + 1 + k)
893 call this%point_neigh(pt_loc_idx)%push(neigh_el)
897 j = j + (2 + num_neigh)
902 deallocate(recv_buffer)
903 call send_buffer%free()
911 type(
mesh_t),
target,
intent(inout) :: this
914 type(
uset_i8_t),
target :: edge_idx, ghost, owner
917 type(mpi_status) :: status
918 type(mpi_request) :: send_req, recv_req
919 integer,
contiguous,
pointer :: p1(:), p2(:), ns_id(:)
920 integer :: i, j, id, ierr, num_edge_glb, edge_offset, num_edge_loc
921 integer :: k, l , shared_offset, glb_nshared, n_glb_id
922 integer(kind=i8) :: C, glb_max, glb_id
923 integer(kind=i8),
pointer :: glb_ptr
924 integer(kind=i8),
allocatable :: recv_buff(:)
925 logical :: shared_edge
927 integer :: max_recv, src, dst, n_recv
931 allocate(this%ddata%local_to_global_edge(this%meds))
933 call edge_idx%init(this%hte%num_entries())
934 call send_buff%init(this%hte%num_entries())
935 call owner%init(this%hte%num_entries())
937 call glb_to_loc%init(32, i)
945 num_edge_glb = 2* this%meds
946 call mpi_allreduce(mpi_in_place, num_edge_glb, 1, &
949 glb_max = int(num_edge_glb,
i8)
951 call non_shared_edges%init(this%hte%num_entries())
953 call it%init(this%hte)
958 k = this%have_point_glb_idx(edge%x(1))
959 l = this%have_point_glb_idx(edge%x(2))
960 p1 => this%point_neigh(k)%array()
961 p2 => this%point_neigh(l)%array()
963 shared_edge = .false.
966 do i = 1, this%point_neigh(k)%size()
967 do j = 1, this%point_neigh(l)%size()
968 if ((p1(i) .eq. p2(j)) .and. &
969 (p1(i) .lt. 0) .and. (p2(j) .lt. 0))
then
979 if (shared_edge)
then
980 glb_id = ((int(edge%x(1),
i8)) + int(edge%x(2),
i8)*c) + glb_max
981 call glb_to_loc%set(glb_id, id)
982 call edge_idx%add(glb_id)
983 call owner%add(glb_id)
984 call send_buff%push(glb_id)
986 call non_shared_edges%push(id)
992 num_edge_loc = non_shared_edges%size()
993 call mpi_exscan(num_edge_loc, edge_offset, 1, &
995 edge_offset = edge_offset + 1
998 ns_id => non_shared_edges%array()
999 do i = 1, non_shared_edges%size()
1001 edge_offset = edge_offset + 1
1009 call mpi_allreduce(send_buff%size(), max_recv, 1, &
1012 call ghost%init(send_buff%size())
1014 allocate(recv_buff(max_recv))
1016 do i = 1,
size(this%neigh_order)
1020 if (this%neigh(src))
then
1021 call mpi_irecv(recv_buff, max_recv, mpi_integer8, &
1025 if (this%neigh(dst))
then
1029 select type(sbarray=>send_buff%data)
1030 type is (
integer(i8))
1031 call mpi_isend(sbarray, send_buff%size(), mpi_integer8, &
1036 if (this%neigh(src))
then
1037 call mpi_wait(recv_req, status, ierr)
1038 call mpi_get_count(status, mpi_integer8, n_recv, ierr)
1041 if ((edge_idx%element(recv_buff(j))) .and. (src .lt.
pe_rank))
then
1042 call ghost%add(recv_buff(j))
1043 call owner%remove(recv_buff(j))
1048 if (this%neigh(dst))
then
1049 call mpi_wait(send_req, mpi_status_ignore, ierr)
1055 glb_nshared = num_edge_loc
1056 call mpi_allreduce(mpi_in_place, glb_nshared, 1, &
1060 call mpi_exscan(owner%size(), shared_offset, 1, &
1062 shared_offset = shared_offset + glb_nshared + 1
1065 call send_buff%clear()
1066 call owner%iter_init()
1067 do while (owner%iter_next())
1068 glb_ptr => owner%iter_value()
1069 if (glb_to_loc%get(glb_ptr, id) .eq. 0)
then
1073 call send_buff%push(glb_ptr)
1074 glb_id = int(shared_offset, i8)
1075 call send_buff%push(glb_id)
1077 shared_offset = shared_offset + 1
1086 this%glb_meds = shared_offset -1
1087 call mpi_allreduce(mpi_in_place, this%glb_meds, 1, &
1094 call mpi_allreduce(send_buff%size(), max_recv, 1, &
1097 deallocate(recv_buff)
1098 allocate(recv_buff(max_recv))
1101 do i = 1,
size(this%neigh_order)
1105 if (this%neigh(src))
then
1106 call mpi_irecv(recv_buff, max_recv, mpi_integer8, &
1110 if (this%neigh(dst))
then
1114 select type(sbarray=>send_buff%data)
1115 type is (
integer(i8))
1116 call mpi_isend(sbarray, send_buff%size(), mpi_integer8, &
1121 if (this%neigh(src))
then
1122 call mpi_wait(recv_req, status, ierr)
1123 call mpi_get_count(status, mpi_integer8, n_recv, ierr)
1126 if (ghost%element(recv_buff(j)))
then
1127 if (glb_to_loc%get(recv_buff(j), id) .eq. 0)
then
1128 n_glb_id = int(recv_buff(j + 1 ), 4)
1138 if (this%neigh(dst))
then
1139 call mpi_wait(send_req, mpi_status_ignore, ierr)
1143 deallocate(recv_buff)
1144 call glb_to_loc%free()
1145 call send_buff%free()
1146 call edge_idx%free()
1147 call non_shared_edges%free()
1155 type(
mesh_t),
target,
intent(inout) :: this
1168 type(mpi_status) :: status
1169 type(mpi_request) :: send_req, recv_req
1170 integer,
allocatable :: recv_buff(:)
1171 integer :: non_shared_facets, shared_facets, facet_offset
1172 integer :: id, glb_nshared, shared_offset, owned_facets
1173 integer :: i, j, ierr, max_recv, src, dst, n_recv
1175 shared_facets = this%ddata%shared_facet%size()
1178 if (this%gdim .eq. 2)
then
1179 allocate(this%ddata%local_to_global_facet(this%meds))
1180 call edge_owner%init(this%meds)
1181 call edge_ghost%init(64, i)
1182 non_shared_facets = this%hte%num_entries() - shared_facets
1184 allocate(this%ddata%local_to_global_facet(this%mfcs))
1185 call face_owner%init(this%mfcs)
1186 call face_ghost%init(64, i)
1187 non_shared_facets = this%htf%num_entries() - shared_facets
1193 call mpi_exscan(non_shared_facets, facet_offset, 1, &
1195 facet_offset = facet_offset + 1
1198 if (this%gdim .eq. 2)
then
1199 call edge_it%init(this%hte)
1200 do while (edge_it%next())
1201 call edge_it%data(id)
1202 edge => edge_it%key()
1203 if (.not. this%ddata%shared_facet%element(id))
then
1206 facet_offset = facet_offset + 1
1208 select type(fmp => this%facet_map)
1210 if (fmp%get(edge, facet_data) .eq. 0)
then
1211 if (facet_data%x(2) .lt. 0)
then
1212 if (abs(facet_data%x(2)) .lt. (this%offset_el + 1))
then
1213 call edge_ghost%set(edge, id)
1215 call edge_owner%push(edge)
1224 owned_facets = edge_owner%size()
1226 call face_it%init(this%htf)
1227 do while (face_it%next())
1228 call face_it%data(id)
1229 face => face_it%key()
1230 if (.not. this%ddata%shared_facet%element(id))
then
1233 facet_offset = facet_offset + 1
1235 select type(fmp => this%facet_map)
1237 if (fmp%get(face, facet_data) .eq. 0)
then
1238 if (facet_data%x(2) .lt. 0)
then
1239 if (abs(facet_data%x(2)) .lt. (this%offset_el + 1))
then
1240 call face_ghost%set(face, id)
1242 call face_owner%push(face)
1251 owned_facets = face_owner%size()
1255 glb_nshared = non_shared_facets
1256 call mpi_allreduce(mpi_in_place, glb_nshared, 1, &
1260 call mpi_exscan(owned_facets, shared_offset, 1, &
1262 shared_offset = shared_offset + glb_nshared + 1
1264 if (this%gdim .eq. 2)
then
1266 if (owned_facets .gt. 32)
then
1267 call send_buff%init(owned_facets)
1269 call send_buff%init()
1272 ed => edge_owner%array()
1273 do i = 1, edge_owner%size()
1274 if (this%hte%get(ed(i), id) .eq. 0)
then
1281 call send_buff%push(ed(i)%x(j))
1283 call send_buff%push(shared_offset)
1285 shared_offset = shared_offset + 1
1291 if (owned_facets .gt. 32)
then
1292 call send_buff%init(owned_facets)
1294 call send_buff%init()
1297 fd => face_owner%array()
1298 do i = 1, face_owner%size()
1299 if (this%htf%get(fd(i), id) .eq. 0)
then
1306 call send_buff%push(fd(i)%x(j))
1308 call send_buff%push(shared_offset)
1310 shared_offset = shared_offset + 1
1319 this%glb_mfcs = shared_offset - 1
1320 call mpi_allreduce(mpi_in_place, this%glb_mfcs, 1, &
1327 call mpi_allreduce(send_buff%size(), max_recv, 1, &
1330 allocate(recv_buff(max_recv))
1333 do i = 1,
size(this%neigh_order)
1337 if (this%neigh(src))
then
1338 call mpi_irecv(recv_buff, max_recv, mpi_integer, &
1342 if (this%neigh(dst))
then
1343 call mpi_isend(send_buff%array(), send_buff%size(), mpi_integer, &
1347 if (this%neigh(src))
then
1348 call mpi_wait(recv_req, status, ierr)
1349 call mpi_get_count(status, mpi_integer, n_recv, ierr)
1351 if (this%gdim .eq. 2)
then
1354 recv_edge = (/recv_buff(j), recv_buff(j+1)/)
1357 if (edge_ghost%get(recv_edge, id) .eq. 0)
then
1365 recv_face = (/recv_buff(j), recv_buff(j+1), &
1366 recv_buff(j+2), recv_buff(j+3) /)
1369 if (face_ghost%get(recv_face, id) .eq. 0)
then
1377 if (this%neigh(dst))
then
1378 call mpi_wait(send_req, mpi_status_ignore, ierr)
1383 if (this%gdim .eq. 2)
then
1384 call edge_owner%free()
1385 call edge_ghost%free()
1387 call face_owner%free()
1388 call face_ghost%free()
1391 call send_buff%free()
1392 deallocate(recv_buff)
1399 class(
mesh_t),
target,
intent(inout) :: this
1400 integer,
value :: el
1401 type(
point_t),
target,
intent(inout) :: p1, p2, p3, p4
1403 integer :: p(4), el_glb_idx
1407 this%lconn = .false.
1410 this%lnumr = .false.
1412 call this%add_point(p1, p(1))
1413 call this%add_point(p2, p(2))
1414 call this%add_point(p3, p(3))
1415 call this%add_point(p4, p(4))
1417 ep => this%elements(el)%e
1418 el_glb_idx = el + this%offset_el
1422 call ep%init(el_glb_idx, &
1423 this%points(p(1)), this%points(p(2)), &
1424 this%points(p(3)), this%points(p(4)))
1435 class(
mesh_t),
target,
intent(inout) :: this
1436 integer,
value :: el
1437 type(
point_t),
target,
intent(inout) :: p1, p2, p3, p4, p5, p6, p7, p8
1439 integer :: p(8), el_glb_idx
1444 this%lconn = .false.
1447 this%lnumr = .false.
1449 call this%add_point(p1, p(1))
1450 call this%add_point(p2, p(2))
1451 call this%add_point(p3, p(3))
1452 call this%add_point(p4, p(4))
1453 call this%add_point(p5, p(5))
1454 call this%add_point(p6, p(6))
1455 call this%add_point(p7, p(7))
1456 call this%add_point(p8, p(8))
1458 ep => this%elements(el)%e
1459 el_glb_idx = el + this%offset_el
1462 call ep%init(el_glb_idx, &
1463 this%points(p(1)), this%points(p(2)), &
1464 this%points(p(3)), this%points(p(4)), &
1465 this%points(p(5)), this%points(p(6)), &
1466 this%points(p(7)), this%points(p(8)))
1475 class(
mesh_t),
intent(inout) :: this
1476 type(
point_t),
intent(inout) :: p
1477 integer,
intent(inout) :: idx
1482 this%max_pts_id =
max(this%max_pts_id, tmp)
1484 if (tmp .le. 0)
then
1488 if (this%htp%get(tmp, idx) .gt. 0)
then
1489 this%mpts = this%mpts + 1
1490 call this%htp%set(tmp, this%mpts)
1491 this%points(this%mpts) = p
1499 class(
mesh_t),
intent(inout) :: this
1503 if (this%htf%get(f, idx) .gt. 0)
then
1504 this%mfcs = this%mfcs + 1
1505 call this%htf%set(f, this%mfcs)
1512 class(
mesh_t),
intent(inout) :: this
1516 if (this%hte%get(e, idx) .gt. 0)
then
1517 this%meds = this%meds + 1
1518 call this%hte%set(e, this%meds)
1525 class(
mesh_t),
intent(inout) :: this
1526 integer,
intent(inout) :: f
1527 integer,
intent(inout) :: e
1529 if (e .gt. this%nelv)
then
1533 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1534 (this%gdim .eq. 3 .and. f .gt. 6))
then
1537 this%facet_type(f, e) = 2
1538 call this%wall%add_facet(f, e)
1544 class(
mesh_t),
intent(inout) :: this
1545 integer,
intent(in) :: e
1546 real(kind=
dp),
dimension(5,12),
intent(in) :: curve_data
1547 integer,
dimension(12),
intent(in) :: curve_type
1549 if (e .gt. this%nelv)
then
1552 if ((this%gdim .eq. 2 .and. sum(curve_type(5:8)) .gt. 0) )
then
1555 call this%curve%add_element(e, curve_data, curve_type)
1562 class(
mesh_t),
intent(inout) :: this
1563 integer,
intent(in) :: f
1564 integer,
intent(in) :: e
1566 if (e .gt. this%nelv)
then
1570 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1571 (this%gdim .eq. 3 .and. f .gt. 6))
then
1574 this%facet_type(f, e) = 2
1575 call this%inlet%add_facet(f, e)
1581 class(
mesh_t),
intent(inout) :: this
1582 integer,
intent(in) :: f
1583 integer,
intent(in) :: e
1584 integer,
intent(in) :: label
1586 if (e .gt. this%nelv)
then
1590 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1591 (this%gdim .eq. 3 .and. f .gt. 6))
then
1594 call this%labeled_zones(label)%add_facet(f, e)
1595 this%facet_type(f,e) = -label
1602 class(
mesh_t),
intent(inout) :: this
1603 integer,
intent(inout) :: f
1604 integer,
intent(inout) :: e
1606 if (e .gt. this%nelv)
then
1610 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1611 (this%gdim .eq. 3 .and. f .gt. 6))
then
1614 this%facet_type(f, e) = 1
1615 call this%outlet_normal%add_facet(f, e)
1622 class(
mesh_t),
intent(inout) :: this
1623 integer,
intent(inout) :: f
1624 integer,
intent(inout) :: e
1626 if (e .gt. this%nelv)
then
1630 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1631 (this%gdim .eq. 3 .and. f .gt. 6))
then
1634 this%facet_type(f, e) = 1
1635 call this%outlet%add_facet(f, e)
1641 class(
mesh_t),
intent(inout) :: this
1642 integer,
intent(inout) :: f
1643 integer,
intent(inout) :: e
1645 if (e .gt. this%nelv)
then
1649 if ((this%gdim .eq. 2 .and. f .gt. 4) .or. &
1650 (this%gdim .eq. 3 .and. f .gt. 6))
then
1653 this%facet_type(f, e) = 2
1654 call this%sympln%add_facet(f, e)
1660 class(
mesh_t),
intent(inout) :: this
1661 integer,
intent(in) :: f
1662 integer,
intent(in) :: e
1663 integer,
intent(in) :: pf
1664 integer,
intent(in) :: pe
1665 integer,
intent(inout) :: pids(4)
1666 integer,
dimension(4) :: org_ids
1668 call this%get_facet_ids(f, e, org_ids)
1669 call this%periodic%add_periodic_facet(f, e, pf, pe, pids, org_ids)
1674 class(
mesh_t),
intent(inout) :: this
1675 integer,
intent(in) :: f
1676 integer,
intent(in) :: e
1677 integer,
intent(inout) :: pids(4)
1682 select type(ele => this%elements(e)%e)
1684 call ele%facet_order(t,f)
1687 call ele%facet_order(t2,f)
1697 class(
mesh_t),
intent(inout) :: this
1703 integer :: org_ids(4), pids(4)
1705 integer,
dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1712 integer,
dimension(2, 4) :: edge_nodes = reshape((/1,3,&
1718 do i = 1, this%periodic%size
1719 e = this%periodic%facet_el(i)%x(2)
1720 f = this%periodic%facet_el(i)%x(1)
1721 pe = this%periodic%p_facet_el(i)%x(2)
1722 pf = this%periodic%p_facet_el(i)%x(1)
1723 pids = this%periodic%p_ids(i)%x
1724 call this%get_facet_ids(f, e, pids)
1725 this%periodic%p_ids(i)%x = pids
1727 do i = 1, this%periodic%size
1728 e = this%periodic%facet_el(i)%x(2)
1729 f = this%periodic%facet_el(i)%x(1)
1730 org_ids = this%periodic%org_ids(i)%x
1731 select type(ele => this%elements(e)%e)
1734 pi => ele%pts(face_nodes(j,f))%p
1735 call pi%set_id(org_ids(j))
1739 pi => ele%pts(edge_nodes(j,f))%p
1740 call pi%set_id(org_ids(j))
1748 class(
mesh_t),
intent(inout) :: this
1749 integer,
intent(in) :: f
1750 integer,
intent(in) :: e
1751 integer,
intent(in) :: pf
1752 integer,
intent(in) :: pe
1753 type(
point_t),
pointer :: pi, pj
1754 real(kind=
dp) :: l(3)
1755 integer :: i, j, id, p_local_idx, match
1758 integer,
dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1765 integer,
dimension(2, 4) :: edge_nodes = reshape((/1,3,&
1771 select type(ele => this%elements(e)%e)
1773 select type(elp => this%elements(pe)%e)
1777 l = l + ele%pts(face_nodes(i,f))%p%x(1:3) - &
1778 elp%pts(face_nodes(i,pf))%p%x(1:3)
1782 pi => ele%pts(face_nodes(i,f))%p
1785 pj => elp%pts(face_nodes(j,pf))%p
1786 if (norm2(pi%x(1:3) - pj%x(1:3) - l) .lt. 1d-7)
then
1787 id = min(pi%id(), pj%id())
1790 p_local_idx = this%get_local(this%points(id))
1794 if ( match .gt. 1)
then
1795 call neko_error(
'Multiple matches when creating periodic ids')
1796 else if (match .eq. 0)
then
1797 call neko_error(
'Cannot find matching periodic point')
1802 select type(elp => this%elements(pe)%e)
1806 l = l + ele%pts(edge_nodes(i,f))%p%x(1:3) - &
1807 elp%pts(edge_nodes(i,pf))%p%x(1:3)
1811 pi => ele%pts(edge_nodes(i,f))%p
1813 pj => elp%pts(edge_nodes(j,pf))%p
1815 if (norm2(pi%x(1:3) - pj%x(1:3) - l) .lt. 1d-7)
then
1816 id = min(pi%id(), pj%id())
1819 p_local_idx = this%get_local(this%points(id))
1830 class(
mesh_t),
intent(inout) :: this
1831 integer,
intent(in) :: f
1832 integer,
intent(in) :: e
1833 integer,
intent(in) :: pf
1834 integer,
intent(in) :: pe
1835 integer,
intent(inout) :: pids(4)
1837 integer :: i, id, p_local_idx
1840 integer,
dimension(4, 6) :: face_nodes = reshape((/1,5,7,3,&
1847 select type(ele => this%elements(e)%e)
1850 pi => ele%pts(face_nodes(i,f))%p
1851 call pi%set_id(pids(i))
1852 call this%add_point(pi, id)
1853 p_local_idx = this%get_local(this%points(id))
1861 class(
mesh_t),
intent(inout) :: this
1862 type(
point_t),
intent(inout) :: p
1869 if (this%htp%get(tmp, local_id) .gt. 0)
then
1870 call neko_error(
'Invalid global id (local point)')
1878 class(
mesh_t),
intent(inout) :: this
1882 if (this%hte%get(e, local_id) .gt. 0)
then
1883 call neko_error(
'Invalid global id (local edge)')
1890 class(
mesh_t),
intent(inout) :: this
1894 if (this%htf%get(f, local_id) .gt. 0)
then
1895 call neko_error(
'Invalid global id (local facet)')
1902 class(
mesh_t),
intent(inout) :: this
1904 integer :: global_id
1906 global_id = this%get_local(e)
1908 if (this%gdim .eq. 2)
then
1910 global_id = this%ddata%local_to_global_facet(global_id)
1914 global_id = this%ddata%local_to_global_edge(global_id)
1922 class(
mesh_t),
intent(inout) :: this
1924 integer :: global_id
1926 global_id = this%get_local_facet(f)
1929 global_id = this%ddata%local_to_global_facet(global_id)
1939 class(
mesh_t),
intent(inout) :: this
1940 integer,
intent(inout) :: index
1943 if (this%htp%get(index, local_id) .eq. 1)
then
1952 class(
mesh_t),
intent(inout) :: this
1953 type(
point_t),
intent(inout) :: p
1954 integer :: local_index
1957 local_index = this%get_local(p)
1958 shared = this%ddata%shared_point%element(local_index)
1966 class(
mesh_t),
intent(inout) :: this
1968 integer :: local_index
1970 local_index = this%get_local(e)
1971 if (this%gdim .eq. 2)
then
1972 shared = this%ddata%shared_facet%element(local_index)
1974 shared = this%ddata%shared_edge%element(local_index)
1980 class(
mesh_t),
intent(inout) :: this
1982 integer :: local_index
1985 local_index = this%get_local(f)
1986 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.
integer, parameter, public log_size
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.
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
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.