build_nodeCoords_for_element Subroutine

  subroutine build_nodeCoords_for_element(mesh,e,nGeo,cornerIDs,flag,edgeCurve,nodeXY)
    !! Fill mesh%nodeCoords(:,:,:,e) for one element. Place the four
    !! corner nodes from the file's node table, then perform linear
    !! (nGeo=1) placement or transfinite interpolation (nGeo>1) using
    !! the (possibly curved) edge curves. For straight sides the edge
    !! curve is filled in here by linear interpolation between corners
    !! at Chebyshev-Gauss-Lobatto parametric coordinates.
    implicit none
    class(Mesh2D_t),intent(inout) :: mesh
    integer,intent(in) :: e
    integer,intent(in) :: nGeo
    integer,intent(in) :: cornerIDs(1:4)
    integer,intent(in) :: flag(1:4)
    real(prec),intent(in) :: edgeCurve(1:2,1:nGeo+1,1:4)
    real(prec),intent(in) :: nodeXY(:,:)
    ! Local
    real(prec) :: P(1:2,1:4)
    real(prec) :: gS(1:2,1:nGeo+1),gE(1:2,1:nGeo+1)
    real(prec) :: gN(1:2,1:nGeo+1),gW(1:2,1:nGeo+1)
    real(prec) :: xi(0:nGeo),wts(0:nGeo)
    real(prec) :: u(1:nGeo+1),t
    integer :: i,j

    P(1:2,1) = nodeXY(1:2,cornerIDs(1))
    P(1:2,2) = nodeXY(1:2,cornerIDs(2))
    P(1:2,3) = nodeXY(1:2,cornerIDs(3))
    P(1:2,4) = nodeXY(1:2,cornerIDs(4))

    if(nGeo == 1) then
      ! Pure bilinear element with the 4 corners
      mesh%nodeCoords(1:2,1,1,e) = P(1:2,1)
      mesh%nodeCoords(1:2,nGeo+1,1,e) = P(1:2,2)
      mesh%nodeCoords(1:2,nGeo+1,nGeo+1,e) = P(1:2,3)
      mesh%nodeCoords(1:2,1,nGeo+1,e) = P(1:2,4)
      return
    endif

    ! Chebyshev-Gauss-Lobatto parametric coordinates on [-1,1] and
    ! their [0,1] image used in the TFI blending weights. We use
    ! SELF's quadrature routine to keep node placement consistent
    ! with the rest of the code.
    call ChebyshevQuadrature(nGeo,xi,wts,CHEBYSHEV_GAUSS_LOBATTO)
    do i = 1,nGeo+1
      u(i) = 0.5_prec*(xi(i-1)+1.0_prec)
    enddo

    ! Edge curves. HOHQMesh's `edgeMap(2,4)` is `(1,2 ; 2,3 ; 4,3 ; 1,4)`,
    ! so its sides 1/2/3/4 already run in SELF's (+ξ, +η, +ξ, +η)
    ! directions. No reversal is needed.
    do i = 1,nGeo+1
      if(flag(1) == 1) then
        gS(1:2,i) = edgeCurve(1:2,i,1)
      else
        t = u(i)
        gS(1:2,i) = (1.0_prec-t)*P(1:2,1)+t*P(1:2,2)
      endif

      if(flag(2) == 1) then
        gE(1:2,i) = edgeCurve(1:2,i,2)
      else
        t = u(i)
        gE(1:2,i) = (1.0_prec-t)*P(1:2,2)+t*P(1:2,3)
      endif

      if(flag(3) == 1) then
        gN(1:2,i) = edgeCurve(1:2,i,3)
      else
        t = u(i)
        gN(1:2,i) = (1.0_prec-t)*P(1:2,4)+t*P(1:2,3)
      endif

      if(flag(4) == 1) then
        gW(1:2,i) = edgeCurve(1:2,i,4)
      else
        t = u(i)
        gW(1:2,i) = (1.0_prec-t)*P(1:2,1)+t*P(1:2,4)
      endif
    enddo

    ! Transfinite (Coons-patch) interpolation. u,v in [0,1].
    do j = 1,nGeo+1
      do i = 1,nGeo+1
        mesh%nodeCoords(1:2,i,j,e) = &
          (1.0_prec-u(j))*gS(1:2,i)+u(j)*gN(1:2,i)+ &
          (1.0_prec-u(i))*gW(1:2,j)+u(i)*gE(1:2,j)- &
          ((1.0_prec-u(i))*(1.0_prec-u(j))*P(1:2,1)+ &
           u(i)*(1.0_prec-u(j))*P(1:2,2)+ &
           u(i)*u(j)*P(1:2,3)+ &
           (1.0_prec-u(i))*u(j)*P(1:2,4))
      enddo
    enddo
  endsubroutine build_nodeCoords_for_element