49 use,
intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr
54 integer(kind=i8),
allocatable :: dof(:,:,:,:)
55 logical,
allocatable :: shared_dof(:,:,:,:)
56 real(kind=
rp),
allocatable :: x(:,:,:,:)
57 real(kind=
rp),
allocatable :: y(:,:,:,:)
58 real(kind=
rp),
allocatable :: z(:,:,:,:)
59 integer,
private :: ntot
67 type(c_ptr) :: x_d = c_null_ptr
68 type(c_ptr) :: y_d = c_null_ptr
69 type(c_ptr) :: z_d = c_null_ptr
83 type(
mesh_t),
target,
intent(inout) :: msh
84 type(
space_t),
target,
intent(inout) :: xh
87 if ((msh%gdim .eq. 3 .and. xh%lz .eq. 1) .or. &
88 (msh%gdim .eq. 2 .and. xh%lz .gt. 1))
then
89 call neko_error(
"Invalid dimension of function space for the given mesh")
98 this%ntot = xh%lx* xh%ly * xh%lz * msh%nelv
104 allocate(this%dof(xh%lx, xh%ly, xh%lz, msh%nelv))
105 allocate(this%shared_dof(xh%lx, xh%ly, xh%lz, msh%nelv))
108 this%shared_dof = .false.
111 if (msh%gdim .eq. 3)
then
124 allocate(this%x(xh%lx, xh%ly, xh%lz, msh%nelv))
125 allocate(this%y(xh%lx, xh%ly, xh%lz, msh%nelv))
126 allocate(this%z(xh%lx, xh%ly, xh%lz, msh%nelv))
152 type(
dofmap_t),
intent(inout) :: this
154 if (
allocated(this%dof))
then
158 if (
allocated(this%shared_dof))
then
159 deallocate(this%shared_dof)
162 if (
allocated(this%x))
then
166 if (
allocated(this%y))
then
170 if (
allocated(this%z))
then
180 if (c_associated(this%x_d))
then
184 if (c_associated(this%y_d))
then
188 if (c_associated(this%z_d))
then
204 integer :: il, jl, ix, iy, iz
205 type(
mesh_t),
pointer :: msh
212 ix = mod(jl - 1, 2) * (xh%lx - 1) + 1
213 iy = (mod(jl - 1, 4)/2) * (xh%ly - 1) + 1
214 iz = ((jl - 1)/4) * (xh%lz - 1) + 1
215 this%dof(ix, iy, iz, il) = int(msh%elements(il)%e%pts(jl)%p%id(),
i8)
216 this%shared_dof(ix, iy, iz, il) = &
217 msh%is_shared(msh%elements(il)%e%pts(jl)%p)
225 type(
mesh_t),
pointer :: msh
230 integer(kind=i8) :: num_dofs_edges(3)
231 integer(kind=i8) :: edge_id, edge_offset
232 logical :: shared_dof
238 num_dofs_edges(1) = int(xh%lx - 2,
i8)
239 num_dofs_edges(2) = int(xh%ly - 2,
i8)
240 num_dofs_edges(3) = int(xh%lz - 2,
i8)
241 edge_offset = int(msh%glb_mpts,
i8) + int(1, 4)
245 select type(ep=>msh%elements(i)%e)
250 call ep%edge_id(edge, 1)
251 shared_dof = msh%is_shared(edge)
252 global_id = msh%get_global(edge)
253 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
255 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i))
then
256 do concurrent(j = 2:xh%lx - 1)
258 this%dof(k, 1, 1, i) = edge_id + (j-2)
259 this%shared_dof(k, 1, 1, i) = shared_dof
262 do concurrent(j = 2:xh%lx - 1)
264 this%dof(k, 1, 1, i) = edge_id + (j-2)
265 this%shared_dof(k, 1, 1, i) = shared_dof
269 call ep%edge_id(edge, 3)
270 shared_dof = msh%is_shared(edge)
271 global_id = msh%get_global(edge)
272 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
273 if(int(edge%x(1),
i8) .ne. this%dof(1,1,xh%lz,i))
then
274 do concurrent(j = 2:xh%lx - 1)
276 this%dof(k, 1, xh%lz, i) = edge_id + (j-2)
277 this%shared_dof(k, 1, xh%lz, i) = shared_dof
280 do concurrent(j = 2:xh%lx - 1)
282 this%dof(k, 1, xh%lz, i) = edge_id + (j-2)
283 this%shared_dof(k, 1, xh%lz, i) = shared_dof
287 call ep%edge_id(edge, 2)
288 shared_dof = msh%is_shared(edge)
289 global_id = msh%get_global(edge)
290 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
291 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,1,i))
then
292 do concurrent(j = 2:xh%lx - 1)
294 this%dof(k, xh%ly, 1, i) = edge_id + (j-2)
295 this%shared_dof(k, xh%ly, 1, i) = shared_dof
298 do concurrent(j = 2:xh%lx - 1)
300 this%dof(k, xh%ly, 1, i) = edge_id + (j-2)
301 this%shared_dof(k, xh%ly, 1, i) = shared_dof
305 call ep%edge_id(edge, 4)
306 shared_dof = msh%is_shared(edge)
307 global_id = msh%get_global(edge)
308 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
309 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,xh%lz,i))
then
310 do concurrent(j = 2:xh%lx - 1)
312 this%dof(k, xh%ly, xh%lz, i) = edge_id + (j-2)
313 this%shared_dof(k, xh%ly, xh%lz, i) = shared_dof
316 do concurrent(j = 2:xh%lx - 1)
318 this%dof(k, xh%ly, xh%lz, i) = edge_id + (j-2)
319 this%shared_dof(k, xh%ly, xh%lz, i) = shared_dof
327 call ep%edge_id(edge, 5)
328 shared_dof = msh%is_shared(edge)
329 global_id = msh%get_global(edge)
330 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
331 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i))
then
332 do concurrent(j = 2:xh%ly - 1)
334 this%dof(1, k, 1, i) = edge_id + (j-2)
335 this%shared_dof(1, k, 1, i) = shared_dof
338 do concurrent(j = 2:xh%ly - 1)
340 this%dof(1, k, 1, i) = edge_id + (j-2)
341 this%shared_dof(1, k, 1, i) = shared_dof
345 call ep%edge_id(edge, 7)
346 shared_dof = msh%is_shared(edge)
347 global_id = msh%get_global(edge)
348 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
349 if(int(edge%x(1),
i8) .ne. this%dof(1,1,xh%lz,i))
then
350 do concurrent(j = 2:xh%ly - 1)
352 this%dof(1, k, xh%lz, i) = edge_id + (j-2)
353 this%shared_dof(1, k, xh%lz, i) = shared_dof
356 do concurrent(j = 2:xh%ly - 1)
358 this%dof(1, k, xh%lz, i) = edge_id + (j-2)
359 this%shared_dof(1, k, xh%lz, i) = shared_dof
363 call ep%edge_id(edge, 6)
364 shared_dof = msh%is_shared(edge)
365 global_id = msh%get_global(edge)
366 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
367 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,1,i))
then
368 do concurrent(j = 2:xh%ly - 1)
370 this%dof(xh%lx, k, 1, i) = edge_id + (j-2)
371 this%shared_dof(xh%lx, k, 1, i) = shared_dof
374 do concurrent(j = 2:xh%ly - 1)
376 this%dof(xh%lx, k, 1, i) = edge_id + (j-2)
377 this%shared_dof(xh%lx, k, 1, i) = shared_dof
381 call ep%edge_id(edge,
i8)
382 shared_dof = msh%is_shared(edge)
383 global_id = msh%get_global(edge)
384 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
385 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,xh%lz,i))
then
386 do concurrent(j = 2:xh%ly - 1)
388 this%dof(xh%lx, k, xh%lz, i) = edge_id + (j-2)
389 this%shared_dof(xh%lx, k, xh%lz, i) = shared_dof
392 do concurrent(j = 2:xh%ly - 1)
394 this%dof(xh%lx, k, xh%lz, i) = edge_id + (j-2)
395 this%shared_dof(xh%lx, k, xh%lz, i) = shared_dof
402 call ep%edge_id(edge, 9)
403 shared_dof = msh%is_shared(edge)
404 global_id = msh%get_global(edge)
405 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
406 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i))
then
407 do concurrent(j = 2:xh%lz - 1)
409 this%dof(1, 1, k, i) = edge_id + (j-2)
410 this%shared_dof(1, 1, k, i) = shared_dof
413 do concurrent(j = 2:xh%lz - 1)
415 this%dof(1, 1, k, i) = edge_id + (j-2)
416 this%shared_dof(1, 1, k, i) = shared_dof
420 call ep%edge_id(edge, 10)
421 shared_dof = msh%is_shared(edge)
422 global_id = msh%get_global(edge)
423 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
424 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,1,i))
then
425 do concurrent(j = 2:xh%lz - 1)
427 this%dof(xh%lx, 1, k, i) = edge_id + (j-2)
428 this%shared_dof(xh%lx, 1, k, i) = shared_dof
431 do concurrent(j = 2:xh%lz - 1)
433 this%dof(xh%lx, 1, k, i) = edge_id + (j-2)
434 this%shared_dof(xh%lx, 1, k, i) = shared_dof
438 call ep%edge_id(edge, 11)
439 shared_dof = msh%is_shared(edge)
440 global_id = msh%get_global(edge)
441 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
442 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,1,i))
then
443 do concurrent(j = 2:xh%lz - 1)
445 this%dof(1, xh%ly, k, i) = edge_id + (j-2)
446 this%shared_dof(1, xh%ly, k, i) = shared_dof
449 do concurrent(j = 2:xh%lz - 1)
451 this%dof(1, xh%ly, k, i) = edge_id + (j-2)
452 this%shared_dof(1, xh%ly, k, i) = shared_dof
456 call ep%edge_id(edge, 12)
457 shared_dof = msh%is_shared(edge)
458 global_id = msh%get_global(edge)
459 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
460 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,xh%ly,1,i))
then
461 do concurrent(j = 2:xh%lz - 1)
463 this%dof(xh%lx, xh%ly, k, i) = edge_id + (j-2)
464 this%shared_dof(xh%lx, xh%ly, k, i) = shared_dof
467 do concurrent(j = 2:xh%lz - 1)
469 this%dof(xh%lx, xh%ly, k, i) = edge_id + (j-2)
470 this%shared_dof(xh%lx, xh%ly, k, i) = shared_dof
477 call ep%facet_id(edge, 3)
478 shared_dof = msh%is_shared(edge)
479 global_id = msh%get_global(edge)
480 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
482 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i))
then
483 do concurrent(j = 2:xh%lx - 1)
485 this%dof(k, 1, 1, i) = edge_id + (j-2)
486 this%shared_dof(k, 1, 1, i) = shared_dof
489 do concurrent(j = 2:xh%lx - 1)
491 this%dof(k, 1, 1, i) = edge_id + (j-2)
492 this%shared_dof(k, 1, 1, i) = shared_dof
496 call ep%facet_id(edge, 4)
497 shared_dof = msh%is_shared(edge)
498 global_id = msh%get_global(edge)
499 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
500 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,1,i))
then
501 do concurrent(j = 2:xh%lx - 1)
503 this%dof(k, xh%ly, 1, i) = edge_id + (j-2)
504 this%shared_dof(k, xh%ly, 1, i) = shared_dof
507 do concurrent(j = 2:xh%lx - 1)
509 this%dof(k, xh%ly, 1, i) = edge_id + (j-2)
510 this%shared_dof(k, xh%ly, 1, i) = shared_dof
517 call ep%facet_id(edge, 1)
518 shared_dof = msh%is_shared(edge)
519 global_id = msh%get_global(edge)
520 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
521 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i))
then
522 do concurrent(j = 2:xh%ly - 1)
524 this%dof(1, k, 1, i) = edge_id + (j-2)
525 this%shared_dof(1, k, 1, i) = shared_dof
528 do concurrent(j = 2:xh%ly - 1)
530 this%dof(1, k, 1, i) = edge_id + (j-2)
531 this%shared_dof(1, k, 1, i) = shared_dof
535 call ep%facet_id(edge, 2)
536 shared_dof = msh%is_shared(edge)
537 global_id = msh%get_global(edge)
538 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
539 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,1,i))
then
540 do concurrent(j = 2:xh%ly - 1)
542 this%dof(xh%lx, k, 1, i) = edge_id + (j-2)
543 this%shared_dof(xh%lx, k, 1, i) = shared_dof
546 do concurrent(j = 2:xh%ly - 1)
548 this%dof(xh%lx, k, 1, i) = edge_id + (j-2)
549 this%shared_dof(xh%lx, k, 1, i) = shared_dof
560 type(
mesh_t),
pointer :: msh
565 integer(kind=i8) :: num_dofs_faces(3)
566 integer(kind=i8) :: facet_offset, facet_id
567 logical :: shared_dof
573 facet_offset = int(msh%glb_mpts,
i8) + &
574 int(msh%glb_meds,
i8) * int(xh%lx-2,
i8) + int(1,
i8)
577 num_dofs_faces(1) = int((xh%ly - 2) * (xh%lz - 2),
i8)
578 num_dofs_faces(2) = int((xh%lx - 2) * (xh%lz - 2),
i8)
579 num_dofs_faces(3) = int((xh%lx - 2) * (xh%ly - 2),
i8)
586 call msh%elements(i)%e%facet_id(face, 1)
587 call msh%elements(i)%e%facet_order(face_order, 1)
588 shared_dof = msh%is_shared(face)
589 global_id = msh%get_global(face)
590 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(1)
591 do concurrent(j = 2:(xh%ly - 1), k = 2:(xh%lz -1))
592 this%dof(1, j, k, i) = &
594 this%shared_dof(1, j, k, i) = shared_dof
597 call msh%elements(i)%e%facet_id(face, 2)
598 call msh%elements(i)%e%facet_order(face_order, 2)
599 shared_dof = msh%is_shared(face)
600 global_id = msh%get_global(face)
601 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(1)
602 do concurrent(j = 2:(xh%ly - 1), k = 2:(xh%lz -1))
603 this%dof(xh%lx, j, k, i) = &
605 this%shared_dof(xh%lx, j, k, i) = shared_dof
612 call msh%elements(i)%e%facet_id(face, 3)
613 call msh%elements(i)%e%facet_order(face_order, 3)
614 shared_dof = msh%is_shared(face)
615 global_id = msh%get_global(face)
616 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(2)
617 do concurrent(j = 2:(xh%lx - 1), k = 2:(xh%lz - 1))
618 this%dof(j, 1, k, i) = &
620 this%shared_dof(j, 1, k, i) = shared_dof
623 call msh%elements(i)%e%facet_id(face, 4)
624 call msh%elements(i)%e%facet_order(face_order, 4)
625 shared_dof = msh%is_shared(face)
626 global_id = msh%get_global(face)
627 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(2)
628 do concurrent(j = 2:(xh%lx - 1), k = 2:(xh%lz - 1))
629 this%dof(j, xh%ly, k, i) = &
631 this%shared_dof(j, xh%ly, k, i) = shared_dof
638 call msh%elements(i)%e%facet_id(face, 5)
639 call msh%elements(i)%e%facet_order(face_order, 5)
640 shared_dof = msh%is_shared(face)
641 global_id = msh%get_global(face)
642 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(3)
643 do concurrent(j = 2:(xh%lx - 1), k = 2:(xh%ly - 1))
644 this%dof(j, k, 1, i) = &
646 this%shared_dof(j, k, 1, i) = shared_dof
649 call msh%elements(i)%e%facet_id(face, 6)
650 call msh%elements(i)%e%facet_order(face_order, 6)
651 shared_dof = msh%is_shared(face)
652 global_id = msh%get_global(face)
653 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(3)
654 do concurrent(j = 2:(xh%lx - 1), k = 2:(xh%ly - 1))
655 this%dof(j, k, xh%lz, i) = &
657 this%shared_dof(j, k, xh%lz, i) = shared_dof
664 pure function dofmap_facetidx(face_order, face, facet_id, k1, j1, lk1, lj1)
result(facet_idx)
666 integer(kind=i8),
intent(in) :: facet_id
667 integer(kind=i8) :: facet_idx
668 integer,
intent(in) :: k1, j1, lk1, lj1
690 if(face_order%x(1) .eq. face%x(1))
then
691 if(face_order%x(2) .lt. face_order%x(4))
then
692 facet_idx = facet_id + j + k*lj
694 facet_idx = facet_id + j*lk + k
696 else if(face_order%x(2) .eq. face%x(1))
then
697 if(face_order%x(3) .lt. face_order%x(1))
then
698 facet_idx = facet_id + lk*(lj-1-j) + k
700 facet_idx = facet_id + (lj-1-j) + k*lj
702 else if(face_order%x(3) .eq. face%x(1))
then
703 if(face_order%x(4) .lt. face_order%x(2))
then
704 facet_idx = facet_id + (lj-1-j) + lj*(lk-1-k)
706 facet_idx = facet_id + lk*(lj-1-j) + (lk-1-k)
708 else if(face_order%x(4) .eq. face%x(1))
then
709 if(face_order%x(1) .lt. face_order%x(3))
then
710 facet_idx = facet_id + lk*j + (lk-1-k)
712 facet_idx = facet_id + j + lj*(lk-1-k)
722 integer :: i, j, el_idx
723 type(
mesh_t),
pointer :: msh
725 real(kind=
rp) :: rp_curve_data(5), curve_data_tot(5,12)
727 integer :: n_edge, curve_type(12)
732 if (msh%gdim .eq. 3)
then
739 call dofmap_xyzlin(xh, msh, msh%elements(i)%e, this%x(1,1,1,i), &
740 this%y(1,1,1,i), this%z(1,1,1,i))
742 do i =1, msh%curve%size
744 el_idx = msh%curve%curve_el(i)%el_idx
745 curve_type = msh%curve%curve_el(i)%curve_type
746 curve_data_tot = msh%curve%curve_el(i)%curve_data
748 if (curve_type(j) .eq. 4)
then
752 if (midpoint .and. xh%lx .gt. 2)
then
754 this%x(1,1,1,el_idx), this%y(1,1,1,el_idx),&
755 this%z(1,1,1,el_idx),curve_type, curve_data_tot)
758 do i =1, msh%curve%size
759 el_idx = msh%curve%curve_el(i)%el_idx
761 if (msh%curve%curve_el(i)%curve_type(j) .eq. 3)
then
762 rp_curve_data = msh%curve%curve_el(i)%curve_data(1:5,j)
764 this%x(1,1,1,el_idx), &
765 this%y(1,1,1,el_idx), &
766 this%z(1,1,1, el_idx), &
767 xh, msh%elements(el_idx)%e, msh%gdim)
771 if (
associated(msh%apply_deform))
then
772 call msh%apply_deform(this%x, this%y, this%z, xh%lx, xh%ly, xh%lz)
785 type(
mesh_t),
pointer,
intent(in) :: msh
786 type(
space_t),
intent(in) :: Xh
788 real(kind=
rp),
intent(inout) :: x(xh%lx, xh%ly, xh%lz), &
789 y(xh%lx, xh%ly, xh%lz), &
790 z(xh%lx, xh%ly, xh%lz)
791 real(kind=
rp) :: xyzb(2,2,2,3), zgml(xh%lx, 3)
792 real(kind=
rp) :: jx(xh%lx*2)
793 real(kind=
rp) :: jxt(xh%lx*2), jyt(xh%lx*2), jzt(xh%lx*2)
794 real(kind=
rp) :: w(4*xh%lx**3), tmp(xh%lx, xh%lx, xh%lx)
795 real(kind=
rp),
dimension(2),
parameter :: zlin = (/-1d0, 1d0/)
803 call copy(zgml(1,1), xh%zg(1,1), xh%lx)
804 call copy(zgml(1,2), xh%zg(1,2), xh%ly)
805 if (msh%gdim .gt. 2)
then
806 call copy(zgml(1,3), xh%zg(1,3), xh%lz)
813 if (msh%gdim .gt. 2)
then
818 call trsp(jx, xh%lx, jxt, 2)
820 if (msh%gdim .eq. 2)
then
824 if (msh%gdim .gt. 2)
then
825 do concurrent(j = 1:msh%gdim)
826 xyzb(1,1,1,j) =
element%pts(1)%p%x(j)
827 xyzb(2,1,1,j) =
element%pts(2)%p%x(j)
828 xyzb(1,2,1,j) =
element%pts(3)%p%x(j)
829 xyzb(2,2,1,j) =
element%pts(4)%p%x(j)
831 xyzb(1,1,2,j) =
element%pts(5)%p%x(j)
832 xyzb(2,1,2,j) =
element%pts(6)%p%x(j)
833 xyzb(1,2,2,j) =
element%pts(7)%p%x(j)
834 xyzb(2,2,2,j) =
element%pts(8)%p%x(j)
837 do concurrent(j = 1:msh%gdim)
838 xyzb(1,1,1,j) =
element%pts(1)%p%x(j)
839 xyzb(2,1,1,j) =
element%pts(2)%p%x(j)
840 xyzb(1,2,1,j) =
element%pts(3)%p%x(j)
841 xyzb(2,2,1,j) =
element%pts(4)%p%x(j)
844 if (msh%gdim .eq. 3)
then
845 call tensr3(tmp, xh%lx, xyzb(1,1,1,1), 2, jx, jyt, jzt, w)
846 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
847 call tensr3(tmp, xh%ly, xyzb(1,1,1,2), 2, jx, jyt, jzt, w)
848 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
849 call tensr3(tmp, xh%lz, xyzb(1,1,1,3), 2, jx, jyt, jzt, w)
850 call copy(z, tmp, xh%lx*xh%ly*xh%lz)
852 call tnsr2d_el(tmp, xh%lx, xyzb(1,1,1,1), 2, jx, jyt)
853 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
854 call tnsr2d_el(tmp, xh%ly, xyzb(1,1,1,2), 2, jx, jyt)
855 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
860 type(
mesh_t),
pointer,
intent(in) :: msh
861 type(
space_t),
intent(in) :: Xh
863 real(kind=
rp),
dimension(Xh%lx, Xh%ly, Xh%lz),
intent(inout) :: x, y, z
864 integer :: curve_type(12), eindx(12)
865 real(kind=
rp) :: curve_data(5,12), x3(3,3,3), y3(3,3,3), z3(3,3,3)
867 real(kind=
rp),
dimension(3),
parameter :: zquad = (/-1d0, 0d0,1d0/)
868 real(kind=
rp) :: zg(3)
869 real(kind=
rp),
dimension(Xh%lx,Xh%lx,Xh%lx) :: tmp
870 real(kind=
rp) :: jx(xh%lx*3)
871 real(kind=
rp) :: jxt(xh%lx*3), jyt(xh%lx*3), jzt(xh%lx*3)
872 real(kind=
rp) :: w(4*xh%lxyz,2)
873 integer :: j, k, n_edges
874 eindx = (/2 , 6 , 8 , 4, &
879 if (msh%gdim .eq. 3)
then
881 call xh3%init(
gll, 3, 3, 3)
884 call xh3%init(
gll, 3, 3)
889 if (curve_type(k) .eq. 4)
then
890 x3(eindx(k),1,1) = curve_data(1,k)
891 y3(eindx(k),1,1) = curve_data(2,k)
892 z3(eindx(k),1,1) = curve_data(3,k)
898 if (msh%gdim .eq. 3)
then
903 call neko_warning(
' m deformation not supported for 2d yet')
911 if (msh%gdim .gt. 2)
then
916 call trsp(jx, xh%lx, jxt, 3)
917 if (msh%gdim .eq. 3)
then
918 call tensr3(tmp, xh%lx, x3, 3, jx, jyt, jzt, w)
919 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
920 call tensr3(tmp, xh%ly, y3, 3, jx, jyt, jzt, w)
921 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
922 call tensr3(tmp, xh%lz, z3, 3, jx, jyt, jzt, w)
923 call copy(z, tmp, xh%lx*xh%ly*xh%lz)
925 call tnsr2d_el(tmp, xh%lx, x3, 3, jx, jyt)
926 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
927 call tnsr2d_el(tmp, xh%ly, y3, 3, jx, jyt)
928 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
940 integer,
intent(in) :: n
941 real(kind=
rp),
intent(inout) :: x(n, n, n)
942 real(kind=
rp),
intent(in) :: zg(n)
943 real(kind=
rp),
intent(inout) :: e(n, n, n)
944 real(kind=
rp),
intent(inout) :: v(n, n, n)
945 integer :: gh_type, ntot, kk, jj, ii, k, j, i
946 real(kind=
rp) :: si, sj, sk, hi, hj, hk
952 do concurrent(i = 1:ntot)
956 do concurrent(i = 1:n, j = 1:n, k = 1:n, ii = 1:n, jj = 1:n:n-1, kk = 1:n:n-1)
957 si = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
958 sj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
959 sk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
960 v(i,j,k) = v(i,j,k) + si*sj*sk*x(ii,jj,kk)
963 if (gh_type .eq. 1)
then
964 do concurrent(i = 1:ntot)
972 do concurrent(i = 1:ntot)
978 do concurrent(i = 1:n, j = 1:n, k = 1:n, jj=1:n:n-1, kk = 1:n:n-1)
979 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
980 hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
981 e(i,j,k) = e(i,j,k) + hj*hk*(x(i,jj,kk)-v(i,jj,kk))
986 do concurrent(i = 1:n, j = 1:n, k = 1:n, ii = 1:n:n-1, kk = 1:n:n-1)
987 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
988 hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
989 e(i,j,k) = e(i,j,k) + hi*hk*(x(ii,j,kk)-v(ii,j,kk))
994 do concurrent(i = 1:n, j = 1:n, k = 1:n, ii = 1:n:n-1, jj = 1:n:n-1)
995 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
996 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
997 e(i,j,k) = e(i,j,k) + hi*hj*(x(ii,jj,k)-v(ii,jj,k))
1000 do concurrent(i = 1:ntot)
1001 e(i,1,1) = e(i,1,1) + v(i,1,1)
1004 if (gh_type .eq. 2)
then
1005 do concurrent(i = 1:ntot)
1013 do concurrent(i = 1:ntot)
1019 do concurrent(i = 1:n, j = 1:n, k = 1:n, ii = 1:n:n-1)
1020 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
1021 v(i,j,k) = v(i,j,k) + hi*(x(ii,j,k)-e(ii,j,k))
1027 do concurrent(i = 1:n, j = 1:n, k = 1:n, jj = 1:n:n-1)
1028 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
1029 v(i,j,k) = v(i,j,k) + hj*(x(i,jj,k)-e(i,jj,k))
1035 do concurrent(i = 1:n, j = 1:n, k = 1:n, kk = 1:n:n-1)
1036 hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
1037 v(i,j,k) = v(i,j,k) + hk*(x(i,j,kk)-e(i,j,kk))
1040 do concurrent(i = 1:ntot)
1041 v(i,1,1) = v(i,1,1) + e(i,1,1)
1051 integer,
intent(in) :: n
1052 real(kind=
rp),
intent(inout) :: x(n, n)
1053 real(kind=
rp),
intent(in) :: zg(n)
1054 real(kind=
rp),
intent(inout) :: e(n, n)
1055 real(kind=
rp),
intent(inout) :: v(n, n)
1056 integer,
intent(in) :: gh_type
1057 integer :: i,j , jj, ii, ntot
1058 real(kind=
rp) :: si, sj, hi, hj
1068 si = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
1069 sj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
1070 v(i,j) = v(i,j) + si*sj*x(ii,jj)
1075 if (gh_type .eq. 1)
then
1088 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
1089 e(i,j) = e(i,j) + hj*(x(i,jj)-v(i,jj))
1099 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
1100 e(i,j) = e(i,j) + hi*(x(ii,j)-v(ii,j))
1105 call add3(x, e, v, ntot)
1112 integer,
intent(in) :: isid, gdim
1113 type(
space_t),
intent(in) :: Xh
1115 real(kind=
rp),
dimension(5),
intent(in) :: curve_data
1116 real(kind=
rp),
dimension(Xh%lx, Xh%ly, Xh%lz),
intent(inout) :: x, y, z
1117 real(kind=
rp) :: pt1x, pt1y, pt2x, pt2y, pt12x, pt12y
1118 real(kind=
rp) :: radius, dtheta, r, xys
1119 real(kind=
rp) :: theta0, xcenn, ycenn, h(xh%lx, 3, 2)
1120 real(kind=
rp) :: xcrved(xh%lx), ycrved(xh%lx), xs, ys
1121 integer :: isid1, ixt, iyt, izt, ix, itmp
1122 integer(i4),
dimension(6),
parameter :: fcyc_to_sym = (/3, 2, 4, 1, 5, 6/)
1123 integer(i4),
dimension(12),
parameter :: ecyc_to_sym = (/1, 6, 2, 5, 3, 8, &
1124 & 4, 7, 9, 10, 12, 11/)
1125 integer,
parameter,
dimension(2, 12) :: edge_nodes = reshape((/1, 2, 3, 4, 5, 6, &
1126 & 7, 8, 1, 3, 2, 4, 5, 7, 6, 8, 1, 5, 2, 6, 3, 7, 4, 8/), (/2,12/))
1133 itmp = ecyc_to_sym(isid)
1136 pt1x =
element%pts(edge_nodes(1,itmp))%p%x(1)
1137 pt1y =
element%pts(edge_nodes(1,itmp))%p%x(2)
1138 pt2x =
element%pts(edge_nodes(2,itmp))%p%x(1)
1139 pt2y =
element%pts(edge_nodes(2,itmp))%p%x(2)
1141 pt1x =
element%pts(edge_nodes(2,itmp))%p%x(1)
1142 pt1y =
element%pts(edge_nodes(2,itmp))%p%x(2)
1143 pt2x =
element%pts(edge_nodes(1,itmp))%p%x(1)
1144 pt2y =
element%pts(edge_nodes(1,itmp))%p%x(2)
1147 radius = curve_data(1)
1151 xys = sqrt(xs**2 + ys**2)
1153 if (abs(2.0 * radius) <= xys * 1.00001) &
1154 &
call neko_error(
'Radius to small for arced element surface')
1156 dtheta = abs(asin(0.5*xys/radius))
1157 pt12x = (pt1x + pt2x)/2.0
1158 pt12y = (pt1y + pt2y)/2.0
1159 xcenn = pt12x - xs/xys * radius*cos(dtheta)
1160 ycenn = pt12y - ys/xys * radius*cos(dtheta)
1161 theta0 = atan2((pt12y-ycenn), (pt12x-xcenn))
1163 isid1 = mod(isid+4-1, 4)+1
1165 if (radius < 0.0) dtheta = -dtheta
1168 if (isid1.gt.2) ixt=xh%lx+1-ix
1170 xcrved(ixt) = xcenn + abs(radius) * cos(theta0 + r*dtheta) &
1171 - ( h(ix,1,1)*pt1x + h(ix,1,2)*pt2x )
1172 ycrved(ixt) = ycenn + abs(radius) * sin(theta0 + r*dtheta) &
1173 - ( h(ix,1,1)*pt1y + h(ix,1,2)*pt2y )
1178 isid1 = fcyc_to_sym(isid1)
1182 if (isid1 .le. 2)
then
1183 call addtnsr(x, h(1,1,ixt), xcrved, h(1,3,izt) &
1184 ,xh%lx, xh%ly, xh%lz)
1185 call addtnsr(y, h(1,1,ixt), ycrved, h(1,3,izt) &
1186 ,xh%lx, xh%ly, xh%lz)
1188 call addtnsr(x, xcrved, h(1,2,iyt), h(1,3,izt) &
1189 ,xh%lx, xh%ly, xh%lz)
1190 call addtnsr(y, ycrved, h(1,2,iyt), h(1,3,izt) &
1191 ,xh%lx, xh%ly, xh%lz)
1196 integer,
intent(in) :: lx, gdim
1197 real(kind=
rp),
intent(inout) :: h(lx, 3, 2)
1198 real(kind=
rp),
intent(in) :: zgml(lx, 3)
1199 integer :: ix, iy, iz
1202 h(ix,1,1) = (1.0_rp - zgml(ix, 1)) * 0.5_rp
1203 h(ix,1,2) = (1.0_rp + zgml(ix, 1)) * 0.5_rp
1207 h(iy,2,1) = (1.0_rp - zgml(iy, 2)) * 0.5_rp
1208 h(iy,2,2) = (1.0_rp + zgml(iy, 2)) * 0.5_rp
1211 if (gdim .eq. 3)
then
1213 h(iz,3,1) = (1.0_rp - zgml(iz, 3)) * 0.5_rp
1214 h(iz,3,2) = (1.0_rp + zgml(iz, 3)) * 0.5_rp
1217 call rone(h(1,3,1), lx)
1218 call rone(h(1,3,2), lx)
Map a Fortran array to a device (allocate and associate)
Copy data between host and device (or device and device)
Device abstraction, common interface for various accelerators.
integer, parameter, public host_to_device
subroutine, public device_free(x_d)
Deallocate memory on the device.
Defines a mapping of the degrees of freedom.
subroutine gh_face_extend_2d(x, zg, n, gh_type, e, v)
Extend 2D faces into interior via gordon hall gh_type: 1 - vertex only 2 - vertex and faces.
subroutine dofmap_xyzquad(Xh, msh, element, x, y, z, curve_type, curve_data)
type(dofmap_t) function dofmap_init(msh, Xh)
subroutine dofmap_generate_xyz(this)
Generate x,y,z-coordinates for all dofs.
subroutine compute_h(h, zgml, gdim, lx)
subroutine arc_surface(isid, curve_data, x, y, z, Xh, element, gdim)
subroutine dofmap_free(this)
Deallocate the dofmap.
subroutine dofmap_number_edges(this)
Assing numbers to dofs on edges.
subroutine dofmap_number_points(this)
Assign numbers to each dofs on points.
subroutine dofmap_xyzlin(Xh, msh, element, x, y, z)
Generate the x, y, z coordinates of the dofs in a signle element, assuming linear element edges.
pure integer function dofmap_size(this)
Return the total number of dofs in the dofmap, lx*ly*lz*nelv.
subroutine dofmap_number_faces(this)
Assign numbers to dofs on faces.
subroutine gh_face_extend_3d(x, zg, n, gh_type, e, v)
Extend faces into interior via gordon hall gh_type: 1 - vertex only 2 - vertex and edges 3 - vertex,...
pure integer(kind=i8) function dofmap_facetidx(face_order, face, facet_id, k1, j1, lk1, lj1)
Get idx for GLL point on face depending on face ordering k and j.
Fast diagonalization methods from NEKTON.
subroutine, public fd_weights_full(xi, x, n, m, c)
Compute finite-difference stencil weights for evaluating derivatives up to order at a point.
Defines a hexahedron element.
subroutine, public rone(a, n)
Set all elements to one.
subroutine, public add3(a, b, c, n)
Vector addition .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
integer, parameter neko_bcknd_device
integer, parameter, public i8
integer, parameter, public i4
integer, parameter, public rp
Global precision used in computations.
Defines a quadrilateral element.
Defines a function space.
integer, parameter, public gll
subroutine, public tensr3(v, nv, u, nu, A, Bt, Ct, w)
Tensor product .
subroutine, public addtnsr(s, h1, h2, h3, nx, ny, nz)
Maps and adds to S a tensor product form of the three functions H1,H2,H3. This is a single element ro...
subroutine, public trsp(a, lda, b, ldb)
Transpose of a rectangular tensor .
subroutine, public tnsr2d_el(v, nv, u, nu, A, Bt)
Computes .
subroutine, public neko_warning(warning_msg)
Base type for an element.
The function space for the SEM solution fields.