anchor.diff:    !----------------------------------------------
     ! find the normal vector of this quadrilateral
     !----------------------------------------------
  -  subroutine findQuadNorm(this,rst)
  +  function findQuadNorm(this)
       class(typeQuad),intent(in)::this
  -    double precision,intent(out)::rst(3)
  +    double precision findQuadNorm(3)
       call find3PNorm(Node(this%NodeInd(1))%Pos(:),Node(this%NodeInd(2))%Pos(:),&
  -    &               Node(this%NodeInd(3))%Pos(:),rst)
  -  end subroutine
  +    &               Node(this%NodeInd(3))%Pos(:),findQuadNorm)
  +  end function
     
     !-------------------------------------
     ! find the center of this tetrahedron
     !-------------------------------------
     ! Note: actually we are finding the average position of the vertices
  -  subroutine findTetPC(this,rst)
  +  function findTetPC(this)
       class(typeTet),intent(in)::this
  -    double precision,intent(out)::rst(3)
  -    rst(:)=(Node(this%NodeInd(1))%Pos(:)+Node(this%NodeInd(2))%Pos(:)&
  -    &      +Node(this%NodeInd(3))%Pos(:)+Node(this%NodeInd(4))%Pos(:))/4d0
  -  end subroutine
  +    double precision findTetPC(3)
  +    findTetPC(:)=(Node(this%NodeInd(1))%Pos(:)+Node(this%NodeInd(2))%Pos(:)&
  +    &            +Node(this%NodeInd(3))%Pos(:)+Node(this%NodeInd(4))%Pos(:))/4d0
  +  end function
     
     !-------------------------------------
     ! find the volume of this tetrahedron
  
positive.code:   ! find the normal vector of this facet
    !--------------------------------------
    ! Note: actually we are finding the average position of the vertices
    subroutine findFacetNorm(this,rst)
      class(typeFacet),intent(in)::this
      double precision,intent(out)::rst(3)
      select case(this%ShapeType)
        case(2)
          call Tri(this%ShapeInd)%findNorm(rst(:))
        case(3)
          call Quad(this%ShapeInd)%findNorm(rst(:))
        case default
          write(*,'(a,i2)'),'ERROR: unknown facet shapeType: ',this%shapeType
          stop
      end select
    end subroutine
    
    !------------------------------------------
    ! get the geometrical entity of this facet
  
negative.code:   !-----------------------------------------
    ! find the normal vector of this triangle
    !-----------------------------------------
    function findTriNorm(this)
      class(typeTri),intent(in)::this
      double precision findTriNorm(3)
      call find3PNorm(Node(this%NodeInd(1))%Pos(:),Node(this%NodeInd(2))%Pos(:),&
      &               Node(this%NodeInd(3))%Pos(:),findTriNorm)
    end function
    
    !---------------------------------------
    ! find the center of this quadrilateral
    !---------------------------------------
    ! Note: actually we are finding the average position of the vertices
    function findQuadPC(this)
      class(typeQuad),intent(in)::this
      double precision findQuadPC(3)
      findQuadPC(:)=(Node(this%NodeInd(1))%Pos(:)+Node(this%NodeInd(2))%Pos(:)&
      &             +Node(this%NodeInd(3))%Pos(:)+Node(this%NodeInd(4))%Pos(:))/4d0
    end function
    
    !-------------------------------------
    ! find the area of this quadrilateral
  