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 center of this facet
    !-------------------------------
    ! Note: actually we are finding the average position of the vertices
    subroutine findFacetPC(this,rst)
      class(typeFacet),intent(in)::this
      double precision,intent(out)::rst(3)
      select case(this%ShapeType)
        case(2)
          call Tri(this%ShapeInd)%findPC(rst(:))
        case(3)
          call Quad(this%ShapeInd)%findPC(rst(:))
        case default
          write(*,'(a,i2)'),'ERROR: unknown facet shapeType: ',this%shapeType
          stop
      end select
    end subroutine
    
    !-----------------------------
    ! find the area of this facet
  
negative.code:   ! find the center of this element
    !---------------------------------
    ! Note: actually we are finding the average position of the vertices
    function findElePC(this)
      class(typeEle),intent(in)::this
      double precision findElePC(3)
      select case(this%ShapeType)
        case(4)
          findElePC(:)=Tet(this%ShapeInd)%findPC()
        case(5)
          findElePC(:)=Hex(this%ShapeInd)%findPC()
        case default
          write(*,'(a,i2)'),'ERROR: unknown element shapeType: ',this%shapeType
          stop
      end select
    end function
    
    !---------------------------------
    ! find the volume of this element
  