subroutine CalculateContravariantBasis_SEMHex(myGeom)
implicit none
class(SEMHex),intent(inout) :: myGeom
! Local
integer :: iEl,i,j,k
real(prec) :: fac
real(prec) :: mag
type(Vector3D) :: xlgradxm,xmgradxl
type(Vector3D) :: curl_xlgradxm,curl_xmgradxl
! Here we use the curl invariant form from Kopriva (2006)
! to calculate the contravariant basis vectors
call xlgradxm%Init(myGeom%x%interp,1,myGeom%x%nElem)
call xmgradxl%Init(myGeom%x%interp,1,myGeom%x%nElem)
call curl_xlgradxm%Init(myGeom%x%interp,1,myGeom%x%nElem)
call curl_xmgradxl%Init(myGeom%x%interp,1,myGeom%x%nElem)
! Ja^{1:3}_1 (n=1, m=2, l=3) First component of the contravariant basis vectors
do iEl = 1,myGeom%nElem
do k = 1,myGeom%dxds%interp%N+1
do j = 1,myGeom%dxds%interp%N+1
do i = 1,myGeom%dxds%interp%N+1
xlgradxm%interior(i,j,k,iel,1,1:3) = myGeom%x%interior(i,j,k,iel,1,3)*myGeom%dxds%interior(i,j,k,iel,1,2,1:3) ! x(...,l)*dxds(...,m,1:3) ; l=3,m=2
xmgradxl%interior(i,j,k,iel,1,1:3) = myGeom%x%interior(i,j,k,iel,1,2)*myGeom%dxds%interior(i,j,k,iel,1,3,1:3) ! x(...,m)*dxds(...,l,1:3) ; l=3,m=2
enddo
enddo
enddo
enddo
call xlgradxm%Curl(curl_xlgradxm%interior)
call xmgradxl%Curl(curl_xmgradxl%interior)
do iEl = 1,myGeom%nElem
do k = 1,myGeom%dxds%interp%N+1
do j = 1,myGeom%dxds%interp%N+1
do i = 1,myGeom%dxds%interp%N+1
! In our convention, dsdx(i,j) is contravariant vector j, component i
! dsdx(...,n,i) = Ja^{i}_{n} = contravariant vector i, component n;
! Here, i = 1:3, and n=1
myGeom%dsdx%interior(i,j,k,iel,1,1,1:3) = 0.5_prec*( &
curl_xmgradxl%interior(i,j,k,iel,1,1:3)- &
curl_xlgradxm%interior(i,j,k,iel,1,1:3))
enddo
enddo
enddo
enddo
! Ja^{1:3}_2 (n=2, m=3, l=1) Second component of the contravariant basis vectors
do iEl = 1,myGeom%nElem
do k = 1,myGeom%dxds%interp%N+1
do j = 1,myGeom%dxds%interp%N+1
do i = 1,myGeom%dxds%interp%N+1
xlgradxm%interior(i,j,k,iel,1,1:3) = myGeom%x%interior(i,j,k,iel,1,1)*myGeom%dxds%interior(i,j,k,iel,1,3,1:3) ! x(...,l)*dxds(...,m,1:3) ; l=1,m=3
xmgradxl%interior(i,j,k,iel,1,1:3) = myGeom%x%interior(i,j,k,iel,1,3)*myGeom%dxds%interior(i,j,k,iel,1,1,1:3) ! x(...,m)*dxds(...,l,1:3) ; l=1,m=3
enddo
enddo
enddo
enddo
call xlgradxm%Curl(curl_xlgradxm%interior)
call xmgradxl%Curl(curl_xmgradxl%interior)
do iEl = 1,myGeom%nElem
do k = 1,myGeom%dxds%interp%N+1
do j = 1,myGeom%dxds%interp%N+1
do i = 1,myGeom%dxds%interp%N+1
! In our convention, dsdx(i,j) is contravariant vector j, component i
! dsdx(...,n,i) = Ja^{i}_{n} = contravariant vector i, component n;
! Here, i = 1:3, and n=2
myGeom%dsdx%interior(i,j,k,iel,1,2,1:3) = 0.5_prec*( &
curl_xmgradxl%interior(i,j,k,iel,1,1:3)- &
curl_xlgradxm%interior(i,j,k,iel,1,1:3))
enddo
enddo
enddo
enddo
! Ja^{1:3}_3 (n=3, m=1, l=2) Third component of the contravariant basis vectors
do iEl = 1,myGeom%nElem
do k = 1,myGeom%dxds%interp%N+1
do j = 1,myGeom%dxds%interp%N+1
do i = 1,myGeom%dxds%interp%N+1
xlgradxm%interior(i,j,k,iel,1,1:3) = myGeom%x%interior(i,j,k,iel,1,2)*myGeom%dxds%interior(i,j,k,iel,1,1,1:3) ! x(...,l)*dxds(...,m,1:3) ; l=2,m=1
xmgradxl%interior(i,j,k,iel,1,1:3) = myGeom%x%interior(i,j,k,iel,1,1)*myGeom%dxds%interior(i,j,k,iel,1,2,1:3) ! x(...,m)*dxds(...,l,1:3) ; l=2,m=1
enddo
enddo
enddo
enddo
call xlgradxm%Curl(curl_xlgradxm%interior)
call xmgradxl%Curl(curl_xmgradxl%interior)
do iEl = 1,myGeom%nElem
do k = 1,myGeom%dxds%interp%N+1
do j = 1,myGeom%dxds%interp%N+1
do i = 1,myGeom%dxds%interp%N+1
! In our convention, dsdx(i,j) is contravariant vector j, component i
! dsdx(...,n,i) = Ja^{i}_{n} = contravariant vector i, component n;
! Here, i = 1:3, and n=3
myGeom%dsdx%interior(i,j,k,iel,1,3,1:3) = 0.5_prec*( &
curl_xmgradxl%interior(i,j,k,iel,1,1:3)- &
curl_xlgradxm%interior(i,j,k,iel,1,1:3))
enddo
enddo
enddo
enddo
call xlgradxm%Free()
call xmgradxl%Free()
call curl_xlgradxm%Free()
call curl_xmgradxl%Free()
! Interpolate the contravariant tensor to the boundaries
call myGeom%dsdx%BoundaryInterp() ! Tensor boundary interp is not offloaded
! Now, calculate nHat (outward pointing normal)
do iEl = 1,myGeom%nElem
do k = 1,6
do j = 1,myGeom%J%interp%N+1
do i = 1,myGeom%J%interp%N+1
if(k == selfSide3D_Top .or. k == selfSide3D_East .or. k == selfSide3D_North) then
fac = sign(1.0_prec,myGeom%J%boundary(i,j,k,iEl,1))
else
fac = -sign(1.0_prec,myGeom%J%boundary(i,j,k,iEl,1))
endif
if(k == 1) then ! Bottom
mag = sqrt(myGeom%dsdx%boundary(i,j,k,iEl,1,1,3)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,2,3)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,3,3)**2)
myGeom%nScale%boundary(i,j,k,iEl,1) = mag
myGeom%nHat%boundary(i,j,k,iEl,1,1:3) = &
fac*myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,3)/mag
! Set the directionality for dsdx on the boundaries
! This is primarily used for DG gradient calculations,
! which do not use nHat for the boundary terms.
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,3) = &
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,3)*fac
elseif(k == 2) then ! South
mag = sqrt(myGeom%dsdx%boundary(i,j,k,iEl,1,1,2)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,2,2)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,3,2)**2)
myGeom%nScale%boundary(i,j,k,iEl,1) = mag
myGeom%nHat%boundary(i,j,k,iEl,1,1:3) = &
fac*myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,2)/mag
! Set the directionality for dsdx on the boundaries
! This is primarily used for DG gradient calculations,
! which do not use nHat for the boundary terms.
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,2) = &
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,2)*fac
elseif(k == 3) then ! East
mag = sqrt(myGeom%dsdx%boundary(i,j,k,iEl,1,1,1)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,2,1)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,3,1)**2)
myGeom%nScale%boundary(i,j,k,iEl,1) = mag
myGeom%nHat%boundary(i,j,k,iEl,1,1:3) = &
fac*myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,1)/mag
! Set the directionality for dsdx on the boundaries
! This is primarily used for DG gradient calculations,
! which do not use nHat for the boundary terms.
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,1) = &
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,1)*fac
elseif(k == 4) then ! North
mag = sqrt(myGeom%dsdx%boundary(i,j,k,iEl,1,1,2)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,2,2)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,3,2)**2)
myGeom%nScale%boundary(i,j,k,iEl,1) = mag
myGeom%nHat%boundary(i,j,k,iEl,1,1:3) = &
fac*myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,2)/mag
! Set the directionality for dsdx on the boundaries
! This is primarily used for DG gradient calculations,
! which do not use nHat for the boundary terms.
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,2) = &
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,2)*fac
elseif(k == 5) then ! West
mag = sqrt(myGeom%dsdx%boundary(i,j,k,iEl,1,1,1)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,2,1)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,3,1)**2)
myGeom%nScale%boundary(i,j,k,iEl,1) = mag
myGeom%nHat%boundary(i,j,k,iEl,1,1:3) = &
fac*myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,1)/mag
! Set the directionality for dsdx on the boundaries
! This is primarily used for DG gradient calculations,
! which do not use nHat for the boundary terms.
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,1) = &
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,1)*fac
elseif(k == 6) then ! Top
mag = sqrt(myGeom%dsdx%boundary(i,j,k,iEl,1,1,3)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,2,3)**2+ &
myGeom%dsdx%boundary(i,j,k,iEl,1,3,3)**2)
myGeom%nScale%boundary(i,j,k,iEl,1) = mag
myGeom%nHat%boundary(i,j,k,iEl,1,1:3) = &
fac*myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,3)/mag
! Set the directionality for dsdx on the boundaries
! This is primarily used for DG gradient calculations,
! which do not use nHat for the boundary terms.
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,3) = &
myGeom%dsdx%boundary(i,j,k,iEl,1,1:3,3)*fac
endif
enddo
enddo
enddo
enddo
call myGeom%dsdx%UpdateDevice()
call myGeom%nHat%UpdateDevice()
call myGeom%nScale%UpdateDevice()
endsubroutine CalculateContravariantBasis_SEMHex