op_obj_ApplyLaplacianDC Subroutine

    impure subroutine op_obj_ApplyLaplacianDC(this,rhs,bcs,varname)
      !> Applies Dirichlet boundary conditions to the RHS of a Laplacian
      ! equation.
      class(op_obj),        intent(inout) :: this                              !! Differential operators utility
      type(eulerian_obj_r), intent(inout) :: rhs                               !! Right hand side
      type(bc_set),         intent(inout) :: bcs                               !! Boundary conditions utility
      character(len=*),     intent(in)    :: varname                           !! Variable name
      ! Work variables
      type(extent_obj)  :: extents
      real(wp), pointer :: val(:,:,:)
      integer           :: n,m,id
      integer           :: i,j,k
      logical           :: found
      integer           :: dir
      integer           :: side


      ! Leave if no regions
      if (.not.allocated(bcs%region)) return

      ! Loop over regions
      do id=1,bcs%count

        ! Check Whether we have a BC for Volume Fraction
        found = bcs%CheckBCExists(bcs%region(id)%name,varname)
        if (.not.found) cycle
        if (bcs%GetBCType(bcs%region(id)%name,varname).ne.BC_DIRICHLET) cycle

        ! Get extents
        extents = bcs%GetExtents(bcs%region(id)%name)

        ! Get direction and side of BC
        call bcs%GetSideDirByRegion(bcs%region(id)%name,side,dir)

        ! Get a pointer to the BC values
        call bcs%GetBCPointer(bcs%region(id)%name,varname,val)

        ! Assuming zero gradient
        associate(lo => extents%lo, hi=>extents%hi)
          if (dir.eq.1) then
            do k=lo(3),hi(3)
              do j=lo(2),hi(2)
                do i=lo(1),hi(1)

                  do n=-this%st+1,this%st
                    do m=-this%st,this%st-1
                      if (side.eq.1 .and. (n+m.gt.0)) then
                        rhs%cell(i,j,k) = rhs%cell(i,j,k)  - this%c_d1dx1 (n,i)*this%c_d1dx1m(m,i+n)*val(i,j,k)
                      end if
                      if (side.eq.0 .and. (n+m.lt.0)) then
                        rhs%cell(i,j,k) = rhs%cell(i,j,k)  - this%c_d1dx1 (n,i)*this%c_d1dx1m(m,i+n)*val(i,j,k)
                      end if
                    end do
                  end do

                end do
              end do
            end do
          end if
          if (dir.eq.2) then
            do k=lo(3),hi(3)
              do j=lo(2),hi(2)
                do i=lo(1),hi(1)

                  do n=-this%st+1,this%st
                    do m=-this%st,this%st-1
                      if (side.eq.1 .and. (n+m.gt.0)) then
                        rhs%cell(i,j,k) = rhs%cell(i,j,k)  - this%c_d1dx2 (n,j)*this%c_d1dx2m(m,j+n)*val(i,j,k)
                      end if
                      if (side.eq.0 .and. (n+m.lt.0)) then
                        rhs%cell(i,j,k) = rhs%cell(i,j,k)  - this%c_d1dx2 (n,j)*this%c_d1dx2m(m,j+n)*val(i,j,k)
                      end if
                    end do
                  end do

                end do
              end do
            end do
          end if
          if (dir.eq.3) then
            do k=lo(3),hi(3)
              do j=lo(2),hi(2)
                do i=lo(1),hi(1)

                  do n=-this%st+1,this%st
                    do m=-this%st,this%st-1
                      if (side.eq.1 .and. (n+m.gt.0)) then
                        rhs%cell(i,j,k) = rhs%cell(i,j,k)  - this%c_d1dx3 (n,k)*this%c_d1dx3m(m,k+n)*val(i,j,k)
                      end if
                      if (side.eq.0 .and. (n+m.lt.0)) then
                        rhs%cell(i,j,k) = rhs%cell(i,j,k)  - this%c_d1dx3 (n,k)*this%c_d1dx3m(m,k+n)*val(i,j,k)
                      end if
                    end do
                  end do

                end do
              end do
            end do
          end if
        end associate

        val => null()
      end do

      return
    end subroutine op_obj_ApplyLaplacianDC