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:         if(n<=grid%nC)then ! internal pairs
            fPF=norm2(grid%pP(:,i)-grid%p(:,n))&
            &   /(norm2(grid%pP(:,i)-grid%p(:,m))+norm2(grid%pP(:,i)-grid%p(:,n)))
            taoF(:,:)=fPF*tao(:,:,m)+(1d0-fPF)*tao(:,:,n)
            flow(:)=grid%aP(i)*matmul(grid%normP(:,i),taoF(:,:))
            dRhou(:,m)=dRhou(:,m)+flow(:)
            dRhou(:,n)=dRhou(:,n)-flow(:)
          else ! boundary pairs
  
negative.code:         end if
        end if
      end do
      deallocate(tau)
    end subroutine
    
  end module
  