NewtonInverse_2D Subroutine

  subroutine NewtonInverse_2D(geometry,iEl,N,xTarget,xi,converged)
    !! Newton inverse-map for the 2D SEM element iEl: solve X(xi) = xTarget for
    !! the reference coordinate xi in R^2, where X is the high-order
    !! interpolant. The Jacobian dX/dxi is taken from geometry%dxds (the
    !! covariant basis tensor), interpolated to xi via Lagrange basis.
    !!
    !! Index convention: dxds%interior(i,j,iEl,1,r,c) = d(x_r)/d(s_c).
    implicit none
    type(SEMQuad),intent(in) :: geometry
    integer,intent(in) :: iEl,N
    real(prec),intent(in) :: xTarget(2)
    real(prec),intent(out) :: xi(2)
    logical,intent(out) :: converged
    ! Local
    integer :: iter,i,j,r,c
    real(prec) :: lS(0:N),lT(0:N)
    real(prec) :: xCur(2),jac(2,2),res(2),delta(2)
    real(prec) :: w,det

    xi(1) = 0.0_prec
    xi(2) = 0.0_prec
    converged = .false.

    do iter = 1,newtonMax
      lS = geometry%x%interp%CalculateLagrangePolynomials(xi(1))
      lT = geometry%x%interp%CalculateLagrangePolynomials(xi(2))

      xCur(1) = 0.0_prec
      xCur(2) = 0.0_prec
      jac(1,1) = 0.0_prec
      jac(1,2) = 0.0_prec
      jac(2,1) = 0.0_prec
      jac(2,2) = 0.0_prec
      do j = 1,N+1
        do i = 1,N+1
          w = lS(i-1)*lT(j-1)
          xCur(1) = xCur(1)+w*geometry%x%interior(i,j,iEl,1,1)
          xCur(2) = xCur(2)+w*geometry%x%interior(i,j,iEl,1,2)
          do c = 1,2
            do r = 1,2
              jac(r,c) = jac(r,c)+w*geometry%dxds%interior(i,j,iEl,1,r,c)
            enddo
          enddo
        enddo
      enddo

      res(1) = xTarget(1)-xCur(1)
      res(2) = xTarget(2)-xCur(2)

      det = jac(1,1)*jac(2,2)-jac(1,2)*jac(2,1)
      if(abs(det) < tiny(1.0_prec)*1.0e3_prec) return ! singular; abandon

      delta(1) = (jac(2,2)*res(1)-jac(1,2)*res(2))/det
      delta(2) = (-jac(2,1)*res(1)+jac(1,1)*res(2))/det
      xi(1) = xi(1)+delta(1)
      xi(2) = xi(2)+delta(2)

      if(max(abs(delta(1)),abs(delta(2))) < newtonTolerance) then
        converged = .true.
        return
      endif

      ! Guard against runaway iterates outside a generous box around [-1,1]^2.
      if(abs(xi(1)) > 5.0_prec .or. abs(xi(2)) > 5.0_prec) return
    enddo

  endsubroutine NewtonInverse_2D