MappedDGDivergence_MappedVector3D Subroutine

  subroutine MappedDGDivergence_MappedVector3D(this,df)
      !! Computes the divergence of a 3-D vector using the weak form
      !! On input, the  attribute of the vector
      !! is assigned and the  attribute is set to the physical
      !! directions of the vector. This method will project the vector
      !! onto the contravariant basis vectors.
    implicit none
    class(MappedVector3D),intent(in) :: this
    type(c_ptr),intent(inout) :: df
    ! Local
    real(prec),pointer :: f_p(:,:,:,:,:,:)
    type(c_ptr) :: fc

    ! Contravariant projection
    call ContravariantProjection_3D_gpu(this%interior_gpu, &
                                        this%geometry%dsdx%interior_gpu,this%interp%N,this%nvar,this%nelem)

    call c_f_pointer(this%interior_gpu,f_p, &
                     [this%interp%N+1,this%interp%N+1,this%interp%N+1,this%nelem,this%nvar,3])

    fc = c_loc(f_p(1,1,1,1,1,1))
    call self_blas_matrixop_dim1_3d(this%interp%dgMatrix_gpu,fc,df, &
                                    this%interp%N,this%interp%N,this%nvar,this%nelem,this%blas_handle)

    fc = c_loc(f_p(1,1,1,1,1,2))
    call self_blas_matrixop_dim2_3d(this%interp%dgMatrix_gpu,fc,df,1.0_c_prec, &
                                    this%interp%N,this%interp%N,this%nvar,this%nelem,this%blas_handle)

    fc = c_loc(f_p(1,1,1,1,1,3))
    call self_blas_matrixop_dim3_3d(this%interp%dgMatrix_gpu,fc,df,1.0_c_prec, &
                                    this%interp%N,this%interp%N,this%nvar,this%nelem,this%blas_handle)
    f_p => null()

    ! Boundary terms
    call DG_BoundaryContribution_3D_gpu(this%interp%bmatrix_gpu,this%interp%qweights_gpu, &
                                        this%boundarynormal_gpu,df,this%interp%N,this%nvar,this%nelem)

    call JacobianWeight_3D_gpu(df,this%geometry%J%interior_gpu,this%interp%N,this%nVar,this%nelem)

  endsubroutine MappedDGDivergence_MappedVector3D