cdifs_obj_PrepareSolverOperatorsVLAP Subroutine

    subroutine cdifs_obj_PrepareSolverOperatorsVLAP(this)
      !> Generates the viscous Laplacian operator and adjusts
      ! at the boundaries.
      use leapBC
      implicit none
      class(cdifs_obj), intent(inout) :: this                                  !! CDIFS solver
      ! Work variables
      integer          :: i,j,k,n,m
      integer          :: stm
      logical          :: found
      type(extent_obj) :: extents
      integer          :: dir
      integer          :: bcdir
      integer          :: side
      integer          :: varstag(3)
      character(len=2) :: varname(3)
      integer          :: ierr
      real(wp),         &
           allocatable :: vGRAD(:,:,:,:,:,:)
      integer          :: shift(3)
      integer          :: m_lo
      integer          :: m_hi

      ! vGRAD = d()_i/dx_j
      ! E.g.: vGRAD(i,j,k,1,2) = dV2/dx1 at i,j,k


      ! Define gradient operators for velocity component in x1-dir,
      ! then build the Laplacian for this velocity component
      associate (lo => this%block%lo, hi => this%block%hi, st => this%op%st, bcs => this%bcs, op =>this%op)

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

        allocate(vGRAD(-st:st,lo(1)-st:hi(1)+st,lo(2)-st:hi(2)+st,lo(3)-st:hi(3)+st,3,3),source = 0.0_wp,stat=ierr)
        if (ierr.ne.0) call this%parallel%Stop('Unable to allocate vGRAD array')


        allocate(this%LAP(-stm:stm,lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),3,3), source = 0.0_wp, stat=ierr)
        if (ierr.ne.0) call this%parallel%Stop('Unable to allocate LAP array')

        ! Define gradient operators for velocity component in x1-dir
        do k=lo(3)-st,hi(3)+st-1
          do j=lo(2)-st,hi(2)+st-1
            do i=lo(1)-st,hi(1)+st-1
              do m=-st+1,st
                ! dV1/dx1 on xm(i),ym(j),zm(k)
                vGRAD(m,i,j,k,1,1) = vGRAD(m,i,j,k,1,1) + this%op%c_d1dx1 (m,i)
                ! dV2/dx2 on xm(i),ym(j),zm(k)
                vGRAD(m,i,j,k,2,2) = vGRAD(m,i,j,k,2,2) + this%op%c_d1dx2 (m,j)
                ! dV3/dx3 on xm(i),ym(j),zm(k)
                vGRAD(m,i,j,k,3,3) = vGRAD(m,i,j,k,3,3) + this%op%c_d1dx3 (m,k)
              end do
            end do
          end do
        end do

        do k=lo(3)-st+1,hi(3)+st
          do j=lo(2)-st+1,hi(2)+st
            do i=lo(1)-st+1,hi(1)+st
              do m=-st,st-1
                ! dV1/dx2 on x(i),y(j),zm(k)
                vGRAD(m,i,j,k,2,1) = vGRAD(m,i,j,k,2,1) + this%op%c_d1dx2m(m,j)
                ! dV1/dx3 on x(i),ym(j),z(k)
                vGRAD(m,i,j,k,3,1) = vGRAD(m,i,j,k,3,1) + this%op%c_d1dx3m(m,k)
                ! dV2/dx1 on x(i),y(j),zm(k)
                vGRAD(m,i,j,k,1,2) = vGRAD(m,i,j,k,1,2) + this%op%c_d1dx1m(m,i)
                ! dV2/dx3 on xm(i),y(j),z(k)
                vGRAD(m,i,j,k,3,2) = vGRAD(m,i,j,k,3,2) + this%op%c_d1dx3m(m,k)
                ! dV3/dx1 on x(i),ym(j),z(k)
                vGRAD(m,i,j,k,1,3) = vGRAD(m,i,j,k,1,3) + this%op%c_d1dx1m(m,i)
                ! dV3/dx2 on x(i),y(j),z(k)
                vGRAD(m,i,j,k,2,3) = vGRAD(m,i,j,k,2,3) + this%op%c_d1dx2m(m,j)
              end do

            end do
          end do
        end do

        ! Boundary conditions
        varname=['V1','V2','V3']
        varstag=[  1,   2,   3]
        do n=1,bcs%count
          ! Get direction and side of BC
          call bcs%GetSideDirByRegion(bcs%region(n)%name,side,bcdir)

          ! Check for velocity component in each direction
          do dir=1,3
            found = bcs%CheckBCExists(bcs%region(n)%name, varname(dir))

            if (found) then
              ! Get extents
              extents = bcs%GetExtents(bcs%region(n)%name)

              ! Check type and apply appropriate BC
              select case ( bcs%GetBCType(bcs%region(n)%name, varname(dir)))
              case (BC_DIRICHLET)

                if (bcdir.eq.varstag(dir)) then

                  select case (side)
                  case (BC_RIGHT)
                    m_lo=1
                    m_hi=st
                  case (BC_LEFT)
                    m_lo=-st
                    m_hi=-1
                  end select

                  shift = 0
                  do m=m_lo,m_hi
                    shift(bcdir) = m
                    do k=extents%lo(3),extents%hi(3)
                      do j=extents%lo(2),extents%hi(2)
                        do i=extents%lo(1),extents%hi(1)
                          vGRAD(:,i+shift(1),j+shift(2),k+shift(3),bcdir,dir) = 0.0_wp
                        end do
                      end do
                    end do
                  end do

                else

                  select case (side)
                  case (BC_RIGHT)
                    m_lo=2
                    m_hi=st
                  case (BC_LEFT)
                    m_lo=-st
                    m_hi=-1
                  end select

                  shift = 0
                  do m=m_lo,m_hi
                    shift(bcdir) = m
                    do k=extents%lo(3),extents%hi(3)
                      do j=extents%lo(2),extents%hi(2)
                        do i=extents%lo(1),extents%hi(1)
                          vGRAD(:,i+shift(1),j+shift(2),k+shift(3),bcdir,dir) = 0.0_wp
                        end do
                      end do
                    end do
                  end do

                end if

              case (BC_NEUMANN)

                if (bcdir.eq.varstag(dir)) then

                  select case (side)
                  case (BC_RIGHT)
                    m_lo=1
                    m_hi=st
                  case (BC_LEFT)
                    m_lo=-st
                    m_hi=-1
                  end select

                  shift = 0
                  do m=m_lo,m_hi
                    shift(bcdir) = m
                    do k=extents%lo(3),extents%hi(3)
                      do j=extents%lo(2),extents%hi(2)
                        do i=extents%lo(1),extents%hi(1)
                          vGRAD(:,i+shift(1),j+shift(2),k+shift(3),bcdir,dir) = 0.0_wp
                        end do
                      end do
                    end do
                  end do

                else

                  select case (side)
                  case (BC_RIGHT)
                    m_lo=1
                    m_hi=st
                  case (BC_LEFT)
                    m_lo=-st
                    m_hi=0
                  end select

                  shift = 0
                  do m=m_lo,m_hi
                    shift(bcdir) = m
                    do k=extents%lo(3),extents%hi(3)
                      do j=extents%lo(2),extents%hi(2)
                        do i=extents%lo(1),extents%hi(1)
                          vGRAD(:,i+shift(1),j+shift(2),k+shift(3),bcdir,dir) = 0.0_wp
                        end do
                      end do
                    end do
                  end do

                end if

              end select
            end if
          end do
        end do

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

              do n=-st+1,st
                do m=-st,st-1

                  ! Laplacian for velocity component in x1-dir
                  this%LAP(n+m,i,j,k,1,1) = this%LAP(n+m,i,j,k,1,1) + op%c_d1dx1m(m,i)*vGRAD(n,i+m,j,k,1,1)
                  this%LAP(n+m,i,j,k,2,1) = this%LAP(n+m,i,j,k,2,1) + op%c_d1dx2 (n,j)*vGRAD(m,i,j+n,k,2,1)
                  this%LAP(n+m,i,j,k,3,1) = this%LAP(n+m,i,j,k,3,1) + op%c_d1dx3 (n,k)*vGRAD(m,i,j,k+n,3,1)

                  ! Laplacian for velocity component in x2-dir
                  this%LAP(n+m,i,j,k,1,2) = this%LAP(n+m,i,j,k,1,2) + op%c_d1dx1 (n,i)*vGRAD(m,i+n,j,k,1,2)
                  this%LAP(n+m,i,j,k,2,2) = this%LAP(n+m,i,j,k,2,2) + op%c_d1dx2m(m,j)*vGRAD(n,i,j+m,k,2,2)
                  this%LAP(n+m,i,j,k,3,2) = this%LAP(n+m,i,j,k,3,2) + op%c_d1dx3 (n,k)*vGRAD(m,i,j,k+n,3,2)

                  ! Laplacian for velocity component in x3-dir
                  this%LAP(n+m,i,j,k,1,3) = this%LAP(n+m,i,j,k,1,3) + op%c_d1dx1 (n,i)*vGRAD(m,i+n,j,k,1,3)
                  this%LAP(n+m,i,j,k,2,3) = this%LAP(n+m,i,j,k,2,3) + op%c_d1dx2 (n,j)*vGRAD(m,i,j+n,k,2,3)
                  this%LAP(n+m,i,j,k,3,3) = this%LAP(n+m,i,j,k,3,3) + op%c_d1dx3m(m,k)*vGRAD(n,i,j,k+m,3,3)

                end do
              end do

            end do
          end do
        end do
      end associate

      ! Clear data
      deallocate(vGRAD)

      return
    end subroutine cdifs_obj_PrepareSolverOperatorsVLAP