CalculateContravariantBasis_SEMQuad Subroutine

  subroutine CalculateContravariantBasis_SEMQuad(myGeom)
    implicit none
    class(SEMQuad),intent(inout) :: myGeom
    ! Local
    integer :: iEl,i,j,k
    real(prec) :: fac
    real(prec) :: mag

    ! Now calculate the contravariant basis vectors
    ! In this convention, dsdx(j,i) is contravariant vector i, component j
    ! To project onto contravariant vector i, dot vector along the first dimension
    do iEl = 1,myGeom%nElem
      do j = 1,myGeom%dxds%interp%N+1
        do i = 1,myGeom%dxds%interp%N+1

          myGeom%dsdx%interior(i,j,iel,1,1,1) = myGeom%dxds%interior(i,j,iel,1,2,2)
          myGeom%dsdx%interior(i,j,iel,1,2,1) = -myGeom%dxds%interior(i,j,iel,1,1,2)
          myGeom%dsdx%interior(i,j,iel,1,1,2) = -myGeom%dxds%interior(i,j,iel,1,2,1)
          myGeom%dsdx%interior(i,j,iel,1,2,2) = myGeom%dxds%interior(i,j,iel,1,1,1)

        enddo
      enddo
    enddo

    ! Interpolate the contravariant tensor to the boundaries
    call myGeom%dsdx%BoundaryInterp() ! Tensor boundary interp is not offloaded

    ! Now, modify the sign of dsdx so that
    ! myGeom % dsdx % boundary is equal to the outward pointing normal vector
    do iEl = 1,myGeom%nElem
      do k = 1,4
        do i = 1,myGeom%J%interp%N+1
          if(k == selfSide2D_East .or. k == selfSide2D_North) then
            fac = sign(1.0_prec,myGeom%J%boundary(i,k,iEl,1))
          else
            fac = -sign(1.0_prec,myGeom%J%boundary(i,k,iEl,1))
          endif

          if(k == 1) then ! South

            mag = sqrt(myGeom%dsdx%boundary(i,k,iEl,1,1,2)**2+ &
                       myGeom%dsdx%boundary(i,k,iEl,1,2,2)**2)

            myGeom%nScale%boundary(i,k,iEl,1) = mag

            myGeom%nHat%boundary(i,k,iEl,1,1:2) = &
              fac*myGeom%dsdx%boundary(i,k,iEl,1,1:2,2)/mag

          elseif(k == 2) then ! East

            mag = sqrt(myGeom%dsdx%boundary(i,k,iEl,1,1,1)**2+ &
                       myGeom%dsdx%boundary(i,k,iEl,1,2,1)**2)

            myGeom%nScale%boundary(i,k,iEl,1) = mag

            myGeom%nHat%boundary(i,k,iEl,1,1:2) = &
              fac*myGeom%dsdx%boundary(i,k,iEl,1,1:2,1)/mag

          elseif(k == 3) then ! North

            mag = sqrt(myGeom%dsdx%boundary(i,k,iEl,1,1,2)**2+ &
                       myGeom%dsdx%boundary(i,k,iEl,1,2,2)**2)

            myGeom%nScale%boundary(i,k,iEl,1) = mag

            myGeom%nHat%boundary(i,k,iEl,1,1:2) = &
              fac*myGeom%dsdx%boundary(i,k,iEl,1,1:2,2)/mag

          elseif(k == 4) then ! West

            mag = sqrt(myGeom%dsdx%boundary(i,k,iEl,1,1,1)**2+ &
                       myGeom%dsdx%boundary(i,k,iEl,1,2,1)**2)

            myGeom%nScale%boundary(i,k,iEl,1) = mag

            myGeom%nHat%boundary(i,k,iEl,1,1:2) = &
              fac*myGeom%dsdx%boundary(i,k,iEl,1,1:2,1)/mag

          endif

          ! Set the directionality for dsdx on the boundaries
          myGeom%dsdx%boundary(i,k,iEl,1,1:2,1:2) = &
            myGeom%dsdx%boundary(i,k,iEl,1,1:2,1:2)*fac

        enddo
      enddo
    enddo

    call myGeom%dsdx%UpdateDevice()
    call myGeom%nHat%UpdateDevice()
    call myGeom%nScale%UpdateDevice()

  endsubroutine CalculateContravariantBasis_SEMQuad