anchor.diff:    ! find the center of this hexahedron
     !------------------------------------
     ! Note: actually we are finding the average position of the vertices
  -  subroutine findHexPC(this,rst)
  +  function findHexPC(this)
       class(typeHex),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(:)&
  -    &      +Node(this%NodeInd(5))%Pos(:)+Node(this%NodeInd(6))%Pos(:)&
  -    &      +Node(this%NodeInd(7))%Pos(:)+Node(this%NodeInd(8))%Pos(:))/8d0
  -  end subroutine
  +    double precision findHexPC(3)
  +    findHexPC(:)=(Node(this%NodeInd(1))%Pos(:)+Node(this%NodeInd(2))%Pos(:)&
  +    &            +Node(this%NodeInd(3))%Pos(:)+Node(this%NodeInd(4))%Pos(:)&
  +    &            +Node(this%NodeInd(5))%Pos(:)+Node(this%NodeInd(6))%Pos(:)&
  +    &            +Node(this%NodeInd(7))%Pos(:)+Node(this%NodeInd(8))%Pos(:))/8d0
  +  end function
     
     !------------------------------------
     ! find the volume of this hexahedron
  
positive.code:   !-----------------------------------------
    ! find the normal vector of this triangle
    !-----------------------------------------
    subroutine findTriNorm(this,rst)
      class(typeTri),intent(in)::this
      double precision,intent(out)::rst(3)
      call find3PNorm(Node(this%NodeInd(1))%Pos(:),Node(this%NodeInd(2))%Pos(:),&
      &               Node(this%NodeInd(3))%Pos(:),rst)
    end subroutine
    
    !---------------------------------------
    ! find the center of this quadrilateral
    !---------------------------------------
    ! Note: actually we are finding the average position of the vertices
    subroutine findQuadPC(this,rst)
      class(typeQuad),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
    
    !-------------------------------------
    ! find the area of this quadrilateral
  
negative.code:   ! find the center of this triangle
    !----------------------------------
    ! Note: actually we are finding the average position of the vertices
    function findTriPC(this)
      class(typeTri),intent(in)::this
      double precision findTriPC(3)
      findTriPC(:)=(Node(this%NodeInd(1))%Pos(:)+Node(this%NodeInd(2))%Pos(:)&
      &            +Node(this%NodeInd(3))%Pos(:))/3d0
    end function
    
    !--------------------------------
    ! find the area of this triangle
  