MappedDivergence_MappedVector3D_t Subroutine

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

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=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,k,iEl,iVar,1)
        Fy = this%interior(ii,j,k,iEl,iVar,2)
        Fz = this%interior(ii,j,k,iEl,iVar,3)
        Fc = this%geometry%dsdx%interior(ii,j,k,iEl,1,1,1)*Fx+ &
             this%geometry%dsdx%interior(ii,j,k,iEl,1,2,1)*Fy+ &
             this%geometry%dsdx%interior(ii,j,k,iEl,1,3,1)*Fz
        dfLoc = dfLoc+this%interp%dMatrix(ii,i)*Fc
      enddo
      dF(i,j,k,iel,ivar) = dfLoc

    enddo

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=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,k,iEl,iVar,1)
        Fy = this%interior(i,ii,k,iEl,iVar,2)
        Fz = this%interior(i,ii,k,iEl,iVar,3)
        Fc = this%geometry%dsdx%interior(i,ii,k,iEl,1,1,2)*Fx+ &
             this%geometry%dsdx%interior(i,ii,k,iEl,1,2,2)*Fy+ &
             this%geometry%dsdx%interior(i,ii,k,iEl,1,3,2)*Fz
        dfLoc = dfLoc+this%interp%dMatrix(ii,j)*Fc
      enddo
      dF(i,j,k,iel,ivar) = (dF(i,j,k,iel,ivar)+dfLoc)

    enddo

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=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,j,ii,iEl,iVar,1)
        Fy = this%interior(i,j,ii,iEl,iVar,2)
        Fz = this%interior(i,j,ii,iEl,iVar,3)
        Fc = this%geometry%dsdx%interior(i,j,ii,iEl,1,1,3)*Fx+ &
             this%geometry%dsdx%interior(i,j,ii,iEl,1,2,3)*Fy+ &
             this%geometry%dsdx%interior(i,j,ii,iEl,1,3,3)*Fz
        dfLoc = dfLoc+this%interp%dMatrix(ii,k)*Fc
      enddo
      dF(i,j,k,iel,ivar) = (dF(i,j,k,iel,ivar)+dfLoc)/this%geometry%J%interior(i,j,k,iEl,1)

    enddo

  endsubroutine MappedDivergence_MappedVector3D_t