op_obj_StrainRate Function

    impure function op_obj_StrainRate(this,U,V,W) result(S)
      !> Computes the strain rate tensor from the velocity field.
      ! Result is on mid points (staggering=0).
      ! Tensor is stored as follows:
      ! S = 0.5*( grad(u) + grad(u)^T )
      !     ( S(1) S(4) S(6) )
      !   = ( S(4) S(2) S(5) )
      !     ( S(6) S(5) S(3) )
      class(op_obj),        intent(in) :: this                                 !! Differential operators utility
      type(Eulerian_obj_r), intent(in) :: U                                    !! Fluid velocity field in 1-dir
      type(Eulerian_obj_r), intent(in) :: V                                    !! Fluid velocity field in 2-dir
      type(Eulerian_obj_r), intent(in) :: W                                    !! Fluid velocity field in 3-dir
      type(eulerian_obj_r)             :: S(6)                                 !! Result
      ! Work variables
      character(str8) :: name
      integer :: i,j,k,n,m,l

      do n=1,6
        call S(n)%Initialize(name,this%block,this%parallel,0)
        S(n) = 0.0_wp
      end do

      associate (lo => this%block%lo,  hi => this%block%hi)
        ! Diagonal elements
        !------------------------------
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              S(1)%cell(i,j,k) = dot_product(this%c_d1dx1(:,i),U%cell(i-this%st+1:i+this%st,j,k))
              S(2)%cell(i,j,k) = dot_product(this%c_d1dx2(:,j),V%cell(i,j-this%st+1:j+this%st,k))
              S(3)%cell(i,j,k) = dot_product(this%c_d1dx3(:,k),W%cell(i,j,k-this%st+1:k+this%st))
            end do
          end do
        end do

        ! Off-diagonal elements
        !------------------------------
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)

             ! Interpolate and differentiate
             ! 0.5*(du/dy+dv/dx)
              do l = -this%st+1,this%st
                do n = -this%st, this%st-1
                  do m = -this%st+1,this%st
                    S(4)%cell(i,j,k) = S(4)%cell(i,j,k)                                                       &
                      + 0.5_wp*this%c_d1dx1(l,i)*this%c_intrp1m(n,i+l)*this%c_intrp2(m,j)*V%cell(i+l+n,j+m,k) &
                      + 0.5_wp*this%c_d1dx2(l,j)*this%c_intrp2m(n,j+l)*this%c_intrp1(m,i)*U%cell(i+m,j+l+n,k)
                  end do
                end do
              end do

              ! 0.5*(dw/dy+dv/dz)
              do l = -this%st+1,this%st
                do n = -this%st, this%st-1
                  do m = -this%st+1,this%st
                    S(5)%cell(i,j,k) = S(5)%cell(i,j,k)                                                       &
                      + 0.5_wp*this%c_d1dx2(l,j)*this%c_intrp2m(n,j+l)*this%c_intrp3(m,k)*W%cell(i,j+l+n,k+m) &
                      + 0.5_wp*this%c_d1dx3(l,k)*this%c_intrp3m(n,k+l)*this%c_intrp2(m,j)*V%cell(i,j+m,k+l+n)
                  end do
                end do
              end do

              ! 0.5*(du/dz+dw/dx)
              do l = -this%st+1,this%st
                do n = -this%st, this%st-1
                  do m = -this%st+1,this%st
                    S(6)%cell(i,j,k) = S(6)%cell(i,j,k)                                                       &
                      + 0.5_wp*this%c_d1dx3(l,k)*this%c_intrp3m(n,k+l)*this%c_intrp1(m,i)*U%cell(i+m,j,k+l+n) &
                      + 0.5_wp*this%c_d1dx1(l,i)*this%c_intrp1m(n,i+l)*this%c_intrp3(m,k)*W%cell(i+l+n,j,k+m)
                  end do
                end do
              end do

            end do
          end do
        end do
      end associate

      do n=1,6
        call S(n)%UpdateGhostCells
      end do

      return
    end function op_obj_StrainRate