MappedDGGradient_MappedScalar2D_t Subroutine

  subroutine MappedDGGradient_MappedScalar2D_t(this,df)
    !!
    implicit none
    class(MappedScalar2D_t),intent(in) :: this
    real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:2)
    ! Local
    integer :: iEl,iVar,i,j,ii,idir
    real(prec) :: dfdx,dfdxb,ja,bfl,bfr

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

      dfdx = 0.0_prec
      do ii = 1,this%N+1
        ja = this%geometry%dsdx%interior(ii,j,iel,1,idir,1)
        dfdx = dfdx+this%interp%dgMatrix(ii,i)*this%interior(ii,j,iel,ivar)*ja
      enddo
      bfl = this%avgboundary(j,4,iel,ivar)*this%geometry%dsdx%boundary(j,4,iel,1,idir,1) ! west
      bfr = this%avgboundary(j,2,iel,ivar)*this%geometry%dsdx%boundary(j,2,iel,1,idir,1) ! east
      dfdxb = (this%interp%bMatrix(i,1)*bfl+this%interp%bMatrix(i,2)*bfr)/this%interp%qweights(i)
      df(i,j,iel,ivar,idir) = dfdx+dfdxb

    enddo

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

      dfdx = 0.0_prec
      do ii = 1,this%N+1
        ja = this%geometry%dsdx%interior(i,ii,iel,1,idir,2)
        dfdx = dfdx+this%interp%dgMatrix(ii,j)*this%interior(i,ii,iel,ivar)*ja
      enddo

      bfl = this%avgboundary(i,1,iel,ivar)*this%geometry%dsdx%boundary(i,1,iel,1,idir,2) ! south
      bfr = this%avgboundary(i,3,iel,ivar)*this%geometry%dsdx%boundary(i,3,iel,1,idir,2) ! north
      dfdxb = (this%interp%bMatrix(j,1)*bfl+this%interp%bMatrix(j,2)*bfr)/this%interp%qweights(j)

      df(i,j,iel,ivar,idir) = (df(i,j,iel,ivar,idir)+dfdx+dfdxb)/this%geometry%J%interior(i,j,iEl,1)

    enddo

  endsubroutine MappedDGGradient_MappedScalar2D_t