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)
257 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i)) k = xh%lx+1-j
258 this%dof(k, 1, 1, i) = edge_id
259 this%shared_dof(k, 1, 1, i) = shared_dof
260 edge_id = edge_id + 1
263 call ep%edge_id(edge, 3)
264 shared_dof = msh%is_shared(edge)
265 global_id = msh%get_global(edge)
266 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
269 if(int(edge%x(1),
i8) .ne. this%dof(1,1,xh%lz,i)) k = xh%lx+1-j
270 this%dof(k, 1, xh%lz, i) = edge_id
271 this%shared_dof(k, 1, xh%lz, i) = shared_dof
272 edge_id = edge_id + 1
275 call ep%edge_id(edge, 2)
276 shared_dof = msh%is_shared(edge)
277 global_id = msh%get_global(edge)
278 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
281 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,1,i)) k = xh%lx+1-j
282 this%dof(k, xh%ly, 1, i) = edge_id
283 this%shared_dof(k, xh%ly, 1, i) = shared_dof
284 edge_id = edge_id + 1
287 call ep%edge_id(edge, 4)
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)
293 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,xh%lz,i)) k = xh%lx+1-j
294 this%dof(k, xh%ly, xh%lz, i) = edge_id
295 this%shared_dof(k, xh%ly, xh%lz, i) = shared_dof
296 edge_id = edge_id + 1
303 call ep%edge_id(edge, 5)
304 shared_dof = msh%is_shared(edge)
305 global_id = msh%get_global(edge)
306 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
309 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i)) k = xh%ly+1-j
310 this%dof(1, k, 1, i) = edge_id
311 this%shared_dof(1, k, 1, i) = shared_dof
312 edge_id = edge_id + 1
315 call ep%edge_id(edge, 7)
316 shared_dof = msh%is_shared(edge)
317 global_id = msh%get_global(edge)
318 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
321 if(int(edge%x(1),
i8) .ne. this%dof(1,1,xh%lz,i)) k = xh%ly+1-j
322 this%dof(1, k, xh%lz, i) = edge_id
323 this%shared_dof(1, k, xh%lz, i) = shared_dof
324 edge_id = edge_id + 1
327 call ep%edge_id(edge, 6)
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)
333 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,1,i)) k = xh%ly+1-j
334 this%dof(xh%lx, k, 1, i) = edge_id
335 this%shared_dof(xh%lx, k, 1, i) = shared_dof
336 edge_id = edge_id + 1
339 call ep%edge_id(edge,
i8)
340 shared_dof = msh%is_shared(edge)
341 global_id = msh%get_global(edge)
342 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
345 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,xh%lz,i)) k = xh%lz+1-j
346 this%dof(xh%lx, k, xh%lz, i) = edge_id
347 this%shared_dof(xh%lx, k, xh%lz, i) = shared_dof
348 edge_id = edge_id + 1
354 call ep%edge_id(edge, 9)
355 shared_dof = msh%is_shared(edge)
356 global_id = msh%get_global(edge)
357 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
360 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i)) k = xh%lz+1-j
361 this%dof(1, 1, k, i) = edge_id
362 this%shared_dof(1, 1, k, i) = shared_dof
363 edge_id = edge_id + 1
366 call ep%edge_id(edge, 10)
367 shared_dof = msh%is_shared(edge)
368 global_id = msh%get_global(edge)
369 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
372 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,1,i)) k = xh%lz+1-j
373 this%dof(xh%lx, 1, k, i) = edge_id
374 this%shared_dof(xh%lx, 1, k, i) = shared_dof
375 edge_id = edge_id + 1
378 call ep%edge_id(edge, 11)
379 shared_dof = msh%is_shared(edge)
380 global_id = msh%get_global(edge)
381 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
384 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,1,i)) k = xh%lz+1-j
385 this%dof(1, xh%ly, k, i) = edge_id
386 this%shared_dof(1, xh%ly, k, i) = shared_dof
387 edge_id = edge_id + 1
390 call ep%edge_id(edge, 12)
391 shared_dof = msh%is_shared(edge)
392 global_id = msh%get_global(edge)
393 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(3)
396 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,xh%ly,1,i)) k = xh%lz+1-j
397 this%dof(xh%lx, xh%ly, k, i) = edge_id
398 this%shared_dof(xh%lx, xh%ly, k, i) = shared_dof
399 edge_id = edge_id + 1
405 call ep%facet_id(edge, 3)
406 shared_dof = msh%is_shared(edge)
407 global_id = msh%get_global(edge)
408 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
412 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i)) k = xh%lx+1-j
413 this%dof(k, 1, 1, i) = edge_id
414 this%shared_dof(k, 1, 1, i) = shared_dof
415 edge_id = edge_id + 1
418 call ep%facet_id(edge, 4)
419 shared_dof = msh%is_shared(edge)
420 global_id = msh%get_global(edge)
421 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(1)
424 if(int(edge%x(1),
i8) .ne. this%dof(1,xh%ly,1,i)) k = xh%lx+1-j
425 this%dof(k, xh%ly, 1, i) = edge_id
426 this%shared_dof(k, xh%ly, 1, i) = shared_dof
427 edge_id = edge_id + 1
433 call ep%facet_id(edge, 1)
434 shared_dof = msh%is_shared(edge)
435 global_id = msh%get_global(edge)
436 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
439 if(int(edge%x(1),
i8) .ne. this%dof(1,1,1,i)) k = xh%ly+1-j
440 this%dof(1, k, 1, i) = edge_id
441 this%shared_dof(1, k, 1, i) = shared_dof
442 edge_id = edge_id + 1
445 call ep%facet_id(edge, 2)
446 shared_dof = msh%is_shared(edge)
447 global_id = msh%get_global(edge)
448 edge_id = edge_offset + int((global_id - 1),
i8) * num_dofs_edges(2)
451 if(int(edge%x(1),
i8) .ne. this%dof(xh%lx,1,1,i)) k = xh%ly+1-j
452 this%dof(xh%lx, k, 1, i) = edge_id
453 this%shared_dof(xh%lx, k, 1, i) = shared_dof
454 edge_id = edge_id + 1
465 type(
mesh_t),
pointer :: msh
470 integer(kind=i8) :: num_dofs_faces(3)
471 integer(kind=i8) :: facet_offset, facet_id
472 logical :: shared_dof
478 facet_offset = int(msh%glb_mpts,
i8) + &
479 int(msh%glb_meds,
i8) * int(xh%lx-2,
i8) + int(1,
i8)
482 num_dofs_faces(1) = int((xh%ly - 2) * (xh%lz - 2),
i8)
483 num_dofs_faces(2) = int((xh%lx - 2) * (xh%lz - 2),
i8)
484 num_dofs_faces(3) = int((xh%lx - 2) * (xh%ly - 2),
i8)
491 call msh%elements(i)%e%facet_id(face, 1)
492 call msh%elements(i)%e%facet_order(face_order, 1)
493 shared_dof = msh%is_shared(face)
494 global_id = msh%get_global(face)
495 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(1)
498 this%dof(1, j, k, i) = &
500 this%shared_dof(1, j, k, i) = shared_dof
504 call msh%elements(i)%e%facet_id(face, 2)
505 call msh%elements(i)%e%facet_order(face_order, 2)
506 shared_dof = msh%is_shared(face)
507 global_id = msh%get_global(face)
508 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(1)
511 this%dof(xh%lx, j, k, i) = &
513 this%shared_dof(xh%lx, j, k, i) = shared_dof
521 call msh%elements(i)%e%facet_id(face, 3)
522 call msh%elements(i)%e%facet_order(face_order, 3)
523 shared_dof = msh%is_shared(face)
524 global_id = msh%get_global(face)
525 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(2)
528 this%dof(j, 1, k, i) = &
530 this%shared_dof(j, 1, k, i) = shared_dof
534 call msh%elements(i)%e%facet_id(face, 4)
535 call msh%elements(i)%e%facet_order(face_order, 4)
536 shared_dof = msh%is_shared(face)
537 global_id = msh%get_global(face)
538 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(2)
541 this%dof(j, xh%ly, k, i) = &
543 this%shared_dof(j, xh%ly, k, i) = shared_dof
551 call msh%elements(i)%e%facet_id(face, 5)
552 call msh%elements(i)%e%facet_order(face_order, 5)
553 shared_dof = msh%is_shared(face)
554 global_id = msh%get_global(face)
555 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(3)
558 this%dof(j, k, 1, i) = &
560 this%shared_dof(j, k, 1, i) = shared_dof
564 call msh%elements(i)%e%facet_id(face, 6)
565 call msh%elements(i)%e%facet_order(face_order, 6)
566 shared_dof = msh%is_shared(face)
567 global_id = msh%get_global(face)
568 facet_id = facet_offset + int((global_id - 1),
i8) * num_dofs_faces(3)
571 this%dof(j, k, xh%lz, i) = &
573 this%shared_dof(j, k, xh%lz, i) = shared_dof
581 function dofmap_facetidx(face_order, face, facet_id, k1, j1, lk1, lj1)
result(facet_idx)
583 integer(kind=i8) :: facet_idx, facet_id
584 integer :: k1, j1, lk1, lj1
606 if(face_order%x(1) .eq. face%x(1))
then
607 if(face_order%x(2) .lt. face_order%x(4))
then
608 facet_idx = facet_id + j + k*lj
610 facet_idx = facet_id + j*lk + k
612 else if(face_order%x(2) .eq. face%x(1))
then
613 if(face_order%x(3) .lt. face_order%x(1))
then
614 facet_idx = facet_id + lk*(lj-1-j) + k
616 facet_idx = facet_id + (lj-1-j) + k*lj
618 else if(face_order%x(3) .eq. face%x(1))
then
619 if(face_order%x(4) .lt. face_order%x(2))
then
620 facet_idx = facet_id + (lj-1-j) + lj*(lk-1-k)
622 facet_idx = facet_id + lk*(lj-1-j) + (lk-1-k)
624 else if(face_order%x(4) .eq. face%x(1))
then
625 if(face_order%x(1) .lt. face_order%x(3))
then
626 facet_idx = facet_id + lk*j + (lk-1-k)
628 facet_idx = facet_id + j + lj*(lk-1-k)
638 integer :: i, j, el_idx
639 type(
mesh_t),
pointer :: msh
641 real(kind=
rp) :: rp_curve_data(5), curve_data_tot(5,12)
643 integer :: n_edge, curve_type(12)
648 if (msh%gdim .eq. 3)
then
655 call dofmap_xyzlin(xh, msh, msh%elements(i)%e, this%x(1,1,1,i), &
656 this%y(1,1,1,i), this%z(1,1,1,i))
658 do i =1, msh%curve%size
660 el_idx = msh%curve%curve_el(i)%el_idx
661 curve_type = msh%curve%curve_el(i)%curve_type
662 curve_data_tot = msh%curve%curve_el(i)%curve_data
664 if (curve_type(j) .eq. 4)
then
668 if (midpoint .and. xh%lx .gt. 2)
then
670 this%x(1,1,1,el_idx), this%y(1,1,1,el_idx),&
671 this%z(1,1,1,el_idx),curve_type, curve_data_tot)
674 do i =1, msh%curve%size
675 el_idx = msh%curve%curve_el(i)%el_idx
677 if (msh%curve%curve_el(i)%curve_type(j) .eq. 3)
then
678 rp_curve_data = msh%curve%curve_el(i)%curve_data(1:5,j)
680 this%x(1,1,1,el_idx), &
681 this%y(1,1,1,el_idx), &
682 this%z(1,1,1, el_idx), &
683 xh, msh%elements(el_idx)%e, msh%gdim)
687 if (
associated(msh%apply_deform))
then
688 call msh%apply_deform(this%x, this%y, this%z, xh%lx, xh%ly, xh%lz)
701 type(
mesh_t),
pointer,
intent(in) :: msh
702 type(
space_t),
intent(in) :: Xh
704 real(kind=
rp),
intent(inout) :: x(xh%lx, xh%ly, xh%lz), &
705 y(xh%lx, xh%ly, xh%lz), &
706 z(xh%lx, xh%ly, xh%lz)
707 real(kind=
rp) :: xyzb(2,2,2,3), zgml(xh%lx, 3)
708 real(kind=
rp) :: jx(xh%lx*2)
709 real(kind=
rp) :: jxt(xh%lx*2), jyt(xh%lx*2), jzt(xh%lx*2)
710 real(kind=
rp) :: w(4*xh%lx**3), tmp(xh%lx, xh%lx, xh%lx)
711 real(kind=
rp),
dimension(2),
parameter :: zlin = (/-1d0, 1d0/)
719 call copy(zgml(1,1), xh%zg(1,1), xh%lx)
720 call copy(zgml(1,2), xh%zg(1,2), xh%ly)
721 if (msh%gdim .gt. 2)
then
722 call copy(zgml(1,3), xh%zg(1,3), xh%lz)
729 if (msh%gdim .gt. 2)
then
734 call trsp(jx, xh%lx, jxt, 2)
736 if (msh%gdim .eq. 2)
then
741 xyzb(1,1,1,j) =
element%pts(1)%p%x(j)
742 xyzb(2,1,1,j) =
element%pts(2)%p%x(j)
743 xyzb(1,2,1,j) =
element%pts(3)%p%x(j)
744 xyzb(2,2,1,j) =
element%pts(4)%p%x(j)
746 if (msh%gdim .gt. 2)
then
747 xyzb(1,1,2,j) =
element%pts(5)%p%x(j)
748 xyzb(2,1,2,j) =
element%pts(6)%p%x(j)
749 xyzb(1,2,2,j) =
element%pts(7)%p%x(j)
750 xyzb(2,2,2,j) =
element%pts(8)%p%x(j)
753 if (msh%gdim .eq. 3)
then
754 call tensr3(tmp, xh%lx, xyzb(1,1,1,1), 2, jx, jyt, jzt, w)
755 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
756 call tensr3(tmp, xh%ly, xyzb(1,1,1,2), 2, jx, jyt, jzt, w)
757 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
758 call tensr3(tmp, xh%lz, xyzb(1,1,1,3), 2, jx, jyt, jzt, w)
759 call copy(z, tmp, xh%lx*xh%ly*xh%lz)
761 call tnsr2d_el(tmp, xh%lx, xyzb(1,1,1,1), 2, jx, jyt)
762 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
763 call tnsr2d_el(tmp, xh%ly, xyzb(1,1,1,2), 2, jx, jyt)
764 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
769 type(
mesh_t),
pointer,
intent(in) :: msh
770 type(
space_t),
intent(in) :: Xh
772 real(kind=
rp),
dimension(Xh%lx, Xh%ly, Xh%lz),
intent(inout) :: x, y, z
773 integer :: curve_type(12), eindx(12)
774 real(kind=
rp) :: curve_data(5,12), x3(3,3,3), y3(3,3,3), z3(3,3,3)
776 real(kind=
rp),
dimension(3),
parameter :: zquad = (/-1d0, 0d0,1d0/)
777 real(kind=
rp) :: zg(3)
778 real(kind=
rp),
dimension(Xh%lx,Xh%lx,Xh%lx) :: tmp
779 real(kind=
rp) :: jx(xh%lx*3)
780 real(kind=
rp) :: jxt(xh%lx*3), jyt(xh%lx*3), jzt(xh%lx*3)
781 real(kind=
rp) :: w(4*xh%lxyz,2)
782 integer :: j, k, n_edges
783 eindx = (/2 , 6 , 8 , 4, &
788 if (msh%gdim .eq. 3)
then
790 call xh3%init(
gll, 3, 3, 3)
793 call xh3%init(
gll, 3, 3)
798 if (curve_type(k) .eq. 4)
then
799 x3(eindx(k),1,1) = curve_data(1,k)
800 y3(eindx(k),1,1) = curve_data(2,k)
801 z3(eindx(k),1,1) = curve_data(3,k)
807 if (msh%gdim .eq. 3)
then
812 call neko_warning(
' m deformation not supported for 2d yet')
820 if (msh%gdim .gt. 2)
then
825 call trsp(jx, xh%lx, jxt, 3)
826 if (msh%gdim .eq. 3)
then
827 call tensr3(tmp, xh%lx, x3, 3, jx, jyt, jzt, w)
828 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
829 call tensr3(tmp, xh%ly, y3, 3, jx, jyt, jzt, w)
830 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
831 call tensr3(tmp, xh%lz, z3, 3, jx, jyt, jzt, w)
832 call copy(z, tmp, xh%lx*xh%ly*xh%lz)
834 call tnsr2d_el(tmp, xh%lx, x3, 3, jx, jyt)
835 call copy(x, tmp, xh%lx*xh%ly*xh%lz)
836 call tnsr2d_el(tmp, xh%ly, y3, 3, jx, jyt)
837 call copy(y, tmp, xh%lx*xh%ly*xh%lz)
849 integer,
intent(in) :: n
850 real(kind=
rp),
intent(inout) :: x(n, n, n)
851 real(kind=
rp),
intent(in) :: zg(n)
852 real(kind=
rp),
intent(inout) :: e(n, n, n)
853 real(kind=
rp),
intent(inout) :: v(n, n, n)
854 integer :: gh_type, ntot, kk, jj, ii, k, j, i
855 real(kind=
rp) :: si, sj, sk, hi, hj, hk
868 si = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
869 sj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
870 sk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
871 v(i,j,k) = v(i,j,k) + si*sj*sk*x(ii,jj,kk)
879 if (gh_type .eq. 1)
then
880 call copy(x, v, ntot)
895 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
896 hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
897 e(i,j,k) = e(i,j,k) + hj*hk*(x(i,jj,kk)-v(i,jj,kk))
911 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
912 hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
913 e(i,j,k) = e(i,j,k) + hi*hk*(x(ii,j,kk)-v(ii,j,kk))
927 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
928 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
929 e(i,j,k) = e(i,j,k) + hi*hj*(x(ii,jj,k)-v(ii,jj,k))
936 call add2(e, v, ntot)
938 if (gh_type .eq. 2)
then
939 call copy(x, e, ntot)
953 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
954 v(i,j,k) = v(i,j,k) + hi*(x(ii,j,k)-e(ii,j,k))
966 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
967 v(i,j,k) = v(i,j,k) + hj*(x(i,jj,k)-e(i,jj,k))
979 hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
980 v(i,j,k) = v(i,j,k) + hk*(x(i,j,kk)-e(i,j,kk))
986 call add2(v, e, ntot)
987 call copy(x, v ,ntot)
995 integer,
intent(in) :: n
996 real(kind=
rp),
intent(inout) :: x(n, n)
997 real(kind=
rp),
intent(in) :: zg(n)
998 real(kind=
rp),
intent(inout) :: e(n, n)
999 real(kind=
rp),
intent(inout) :: v(n, n)
1000 integer,
intent(in) :: gh_type
1001 integer :: i,j , jj, ii, ntot
1002 real(kind=
rp) :: si, sj, hi, hj
1012 si = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
1013 sj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
1014 v(i,j) = v(i,j) + si*sj*x(ii,jj)
1019 if (gh_type .eq. 1)
then
1032 hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
1033 e(i,j) = e(i,j) + hj*(x(i,jj)-v(i,jj))
1043 hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
1044 e(i,j) = e(i,j) + hi*(x(ii,j)-v(ii,j))
1049 call add3(x, e, v, ntot)
1056 integer,
intent(in) :: isid, gdim
1057 type(
space_t),
intent(in) :: Xh
1059 real(kind=
rp),
dimension(5),
intent(in) :: curve_data
1060 real(kind=
rp),
dimension(Xh%lx, Xh%ly, Xh%lz),
intent(inout) :: x, y, z
1061 real(kind=
rp) :: pt1x, pt1y, pt2x, pt2y, pt12x, pt12y
1062 real(kind=
rp) :: radius, dtheta, r, xys
1063 real(kind=
rp) :: theta0, xcenn, ycenn, h(xh%lx, 3, 2)
1064 real(kind=
rp) :: xcrved(xh%lx), ycrved(xh%lx), xs, ys
1065 integer :: isid1, ixt, iyt, izt, ix, itmp
1066 integer(i4),
dimension(6),
parameter :: fcyc_to_sym = (/3, 2, 4, 1, 5, 6/)
1067 integer(i4),
dimension(12),
parameter :: ecyc_to_sym = (/1, 6, 2, 5, 3, 8, &
1068 & 4, 7, 9, 10, 12, 11/)
1069 integer,
parameter,
dimension(2, 12) :: edge_nodes = reshape((/1, 2, 3, 4, 5, 6, &
1070 & 7, 8, 1, 3, 2, 4, 5, 7, 6, 8, 1, 5, 2, 6, 3, 7, 4, 8/), (/2,12/))
1077 itmp = ecyc_to_sym(isid)
1080 pt1x =
element%pts(edge_nodes(1,itmp))%p%x(1)
1081 pt1y =
element%pts(edge_nodes(1,itmp))%p%x(2)
1082 pt2x =
element%pts(edge_nodes(2,itmp))%p%x(1)
1083 pt2y =
element%pts(edge_nodes(2,itmp))%p%x(2)
1085 pt1x =
element%pts(edge_nodes(2,itmp))%p%x(1)
1086 pt1y =
element%pts(edge_nodes(2,itmp))%p%x(2)
1087 pt2x =
element%pts(edge_nodes(1,itmp))%p%x(1)
1088 pt2y =
element%pts(edge_nodes(1,itmp))%p%x(2)
1091 radius = curve_data(1)
1095 xys = sqrt(xs**2 + ys**2)
1097 if (abs(2.0 * radius) <= xys * 1.00001) &
1098 &
call neko_error(
'Radius to small for arced element surface')
1100 dtheta = abs(asin(0.5*xys/radius))
1101 pt12x = (pt1x + pt2x)/2.0
1102 pt12y = (pt1y + pt2y)/2.0
1103 xcenn = pt12x - xs/xys * radius*cos(dtheta)
1104 ycenn = pt12y - ys/xys * radius*cos(dtheta)
1105 theta0 = atan2((pt12y-ycenn), (pt12x-xcenn))
1107 isid1 = mod(isid+4-1, 4)+1
1109 if (radius < 0.0) dtheta = -dtheta
1112 if (isid1.gt.2) ixt=xh%lx+1-ix
1114 xcrved(ixt) = xcenn + abs(radius) * cos(theta0 + r*dtheta) &
1115 - ( h(ix,1,1)*pt1x + h(ix,1,2)*pt2x )
1116 ycrved(ixt) = ycenn + abs(radius) * sin(theta0 + r*dtheta) &
1117 - ( h(ix,1,1)*pt1y + h(ix,1,2)*pt2y )
1122 isid1 = fcyc_to_sym(isid1)
1126 if (isid1 .le. 2)
then
1127 call addtnsr(x, h(1,1,ixt), xcrved, h(1,3,izt) &
1128 ,xh%lx, xh%ly, xh%lz)
1129 call addtnsr(y, h(1,1,ixt), ycrved, h(1,3,izt) &
1130 ,xh%lx, xh%ly, xh%lz)
1132 call addtnsr(x, xcrved, h(1,2,iyt), h(1,3,izt) &
1133 ,xh%lx, xh%ly, xh%lz)
1134 call addtnsr(y, ycrved, h(1,2,iyt), h(1,3,izt) &
1135 ,xh%lx, xh%ly, xh%lz)
1140 integer,
intent(in) :: lx, gdim
1141 real(kind=
rp),
intent(inout) :: h(lx, 3, 2)
1142 real(kind=
rp),
intent(in) :: zgml(lx, 3)
1143 integer :: ix, iy, iz
1146 h(ix,1,1) = (1.0_rp - zgml(ix, 1)) * 0.5_rp
1147 h(ix,1,2) = (1.0_rp + zgml(ix, 1)) * 0.5_rp
1151 h(iy,2,1) = (1.0_rp - zgml(iy, 2)) * 0.5_rp
1152 h(iy,2,2) = (1.0_rp + zgml(iy, 2)) * 0.5_rp
1155 if (gdim .eq. 3)
then
1157 h(iz,3,1) = (1.0_rp - zgml(iz, 3)) * 0.5_rp
1158 h(iz,3,2) = (1.0_rp + zgml(iz, 3)) * 0.5_rp
1161 call rone(h(1,3,1), lx)
1162 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.
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.
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,...
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 add2(a, b, 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 neko_warning(warning_msg)
Base type for an element.
The function space for the SEM solution fields.