anchor.diff:          rf(2)=-grid%normP(1,i)*tf(3)+grid%normP(3,i)*tf(1)
           rf(3)=grid%normP(1,i)*tf(2)-grid%normP(2,i)*tf(1)
           Afs=grid%aP(i)/dot_product(sf,grid%normP(:,i))
  -        flow(:)=dF(:)*Afs*((s(:,n)-s(:,m))/dPF&
  -        &                  -matmul(transpose(gradSF(:,:)),&
  -        &                          dot_product(sf,tf)*tf+dot_product(sf,rf)*rf))
           if(m<=grid%nC.and.n<=grid%nC)then ! internal pairs
  +          flow(:)=dF(:)*Afs*((s(:,n)-s(:,m))/dPF&
  +          &                  -matmul(transpose(gradSF(:,:)),&
  +          &                          dot_product(sf,tf)*tf+dot_product(sf,rf)*rf))
             diff(:,m)=diff(:,m)+flow(:)
             diff(:,n)=diff(:,n)-flow(:)
  -        else
  +        else ! boundary pairs
  +          flow(:)=dF(:)*Afs*(s(:,n)-s(:,m))/(2d0*dPF)
             diff(:,m)=diff(:,m)+flow(:)
           end if
         end if
  
positive.code:     deallocate(diffv)
    end subroutine
    
  end module
  
negative.code:     module procedure::findDiffPolyScal
    end interface
    public::findDiff
    
    !> generic finding viscose force
    interface findViscForce
      module procedure::findViscForcePoly
    end interface
    public::findViscForce
    
  contains
    
  