cdifs_obj_PrepareSolverOperatorsPLAP Subroutine

    subroutine cdifs_obj_PrepareSolverOperatorsPLAP(this)
      !> Generates the pressure Laplacian operator.
      implicit none
      class(cdifs_obj), intent(inout) :: this                                  !! CDIFS solver
      ! Work variables
      integer          :: i,j,k,n,m
      integer          :: stm
      integer          :: ierr

      ! Use differential operators to build the laplacian inside HYPRE
      if (allocated(this%hypre%mat)) deallocate(this%hypre%mat)

      associate (lo => this%block%lo, hi => this%block%hi, op => this%op)

        ! Set stencil extent based on current scheme order
        stm = 2*op%st-1
        this%hypre%stm = stm

        allocate(this%hypre%mat(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),3,-stm:stm),&
                 source = 0.0_wp, stat=ierr)
        if (ierr.ne.0) call this%parallel%Stop('Unable to allocate HYPRE%MAT array')
      end associate


      ! Build Laplacian by convolving divergence and gradient operators
      associate (lo => this%block%lo, hi => this%block%hi, op => this%op)

        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)

              do n=-op%st+1,op%st
                do m=-op%st,op%st-1
                  this%hypre%mat(i,j,k,1,n+m) = this%hypre%mat(i,j,k,1,n+m) + this%DIV(n,i,j,k,1)*this%pGRAD(m,i+n,j,k,1)
                  this%hypre%mat(i,j,k,2,n+m) = this%hypre%mat(i,j,k,2,n+m) + this%DIV(n,i,j,k,2)*this%pGRAD(m,i,j+n,k,2)
                  this%hypre%mat(i,j,k,3,n+m) = this%hypre%mat(i,j,k,3,n+m) + this%DIV(n,i,j,k,3)*this%pGRAD(m,i,j,k+n,3)
                end do
              end do

            end do
          end do
        end do

      end associate

      ! Setup HYPRE system of equations
      call this%hypre%setup()

      return
    end subroutine cdifs_obj_PrepareSolverOperatorsPLAP