anchor.diff:    ! find the center of this element
     !---------------------------------
     ! Note: actually we are finding the average position of the vertices
  -  subroutine findElePC(this,rst)
  +  function findElePC(this)
       class(typeEle),intent(in)::this
  -    double precision,intent(out)::rst(3)
  +    double precision findElePC(3)
       select case(this%ShapeType)
         case(4)
  -        call Tet(this%ShapeInd)%findPC(rst(:))
  +        findElePC(:)=Tet(this%ShapeInd)%findPC()
         case(5)
  -        call Hex(this%ShapeInd)%findPC(rst(:))
  +        findElePC(:)=Hex(this%ShapeInd)%findPC()
         case default
           write(*,'(a,i2)'),'ERROR: unknown element shapeType: ',this%shapeType
           stop
       end select
  -  end subroutine
  +  end function
     
     !---------------------------------
     ! find the volume of this element
  
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 quadrilateral
    !----------------------------------------------
    function findQuadNorm(this)
      class(typeQuad),intent(in)::this
      double precision findQuadNorm(3)
      call find3PNorm(Node(this%NodeInd(1))%Pos(:),Node(this%NodeInd(2))%Pos(:),&
      &               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
    function findTetPC(this)
      class(typeTet),intent(in)::this
      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
  