anchor.diff:    end subroutine
     
     !> find Newtonian viscous force on polyFvGrid
  -  !> \f[ \int_A \hat{n} \cdot \mathbf{\tao} dA \f]
  +  !> \f[ \int_A \hat{n} \cdot \mathbf{\tau} dA \f]
     subroutine findViscForcePoly(grid,u,visc,dRhou)
       use modPolyFvGrid
       use modGradient
  
positive.code:     double precision,intent(in)::u(:,:) !< state variables
      double precision,intent(in)::visc(:) !< viscosity
      double precision,allocatable,intent(inout)::dRhou(:,:) !< force (net flow of rhou)
      double precision::fPF,taoF(DIMS,DIMS),flow(DIMS),utP(DIMS),utF(DIMS)
      double precision,allocatable::tao(:,:,:)
      
      call grid%up()
      if(.not.(allocated(dRhou)))then
        allocate(dRhou(DIMS,grid%nC))
      end if
      dRhou(:,:)=0d0
      allocate(tao(DIMS,DIMS,grid%nC))
      call findGrad(grid,u(:,1:grid%nC),tao)
      forall(i=1:grid%nC)
        tao(:,:,i)=visc(i)*((tao(:,:,i)+transpose(tao(:,:,i)))&
        &                   -2d0/3d0*reshape([1d0,0d0,0d0,0d0,1d0,0d0,0d0,0d0,1d0],[DIMS,DIMS]))
      end forall
      do i=1,grid%nP
        m=grid%iEP(1,i)
  
negative.code:         end if
        end if
      end do
      deallocate(tau)
    end subroutine
    
  end module
  