MappedDivergence_MappedVector2D_t Subroutine

  subroutine MappedDivergence_MappedVector2D_t(this,df)
    ! Strong Form Operator
    !    !
    implicit none
    class(MappedVector2D_t),intent(in) :: this
    real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar)
    ! Local
    integer :: iEl,iVar,i,j,ii
    real(prec) :: dfLoc,Fx,Fy,Fc

    do concurrent(i=1:this%N+1,j=1:this%N+1,iel=1:this%nElem,ivar=1:this%nVar)

      dfLoc = 0.0_prec
      do ii = 1,this%N+1
        ! Convert from physical to computational space
        Fx = this%interior(ii,j,iEl,iVar,1)
        Fy = this%interior(ii,j,iEl,iVar,2)
        Fc = this%geometry%dsdx%interior(ii,j,iEl,1,1,1)*Fx+ &
             this%geometry%dsdx%interior(ii,j,iEl,1,2,1)*Fy
        dfLoc = dfLoc+this%interp%dMatrix(ii,i)*Fc
      enddo
      dF(i,j,iel,ivar) = dfLoc

    enddo

    do concurrent(i=1:this%N+1,j=1:this%N+1,iel=1:this%nElem,ivar=1:this%nVar)

      dfLoc = 0.0_prec
      do ii = 1,this%N+1
        ! Convert from physical to computational space
        Fx = this%interior(i,ii,iEl,iVar,1)
        Fy = this%interior(i,ii,iEl,iVar,2)
        Fc = this%geometry%dsdx%interior(i,ii,iEl,1,1,2)*Fx+ &
             this%geometry%dsdx%interior(i,ii,iEl,1,2,2)*Fy
        dfLoc = dfLoc+this%interp%dMatrix(ii,j)*Fc
      enddo
      dF(i,j,iel,ivar) = (dF(i,j,iel,ivar)+dfLoc)/this%geometry%J%interior(i,j,iEl,1)

    enddo

  endsubroutine MappedDivergence_MappedVector2D_t