anchor.diff:            maskHex(:)=.false.
             lPHex(:)=Hex(Ele(i)%ShapeInd)%NodeInd(:)
             do j=1,8
  -            if(this%NodeInd(j)==lPHex(1).or.this%NodeInd(j)==lPHex(2).or.&
  -            &  this%NodeInd(j)==lPHex(3).or.this%NodeInd(j)==lPHex(4).or.&
  -            &  this%NodeInd(j)==lPHex(5).or.this%NodeInd(j)==lPHex(6).or.&
  -            &  this%NodeInd(j)==lPHex(7).or.this%NodeInd(j)==lPHex(8))then
  +            if(any(lPHex(:)==this%NodeInd(j)))then
                 maskHex(j)=.true.
               end if
             end do
  
positive.code:           maskTet(:)=.false.
            lPTet(:)=Tet(Ele(i)%ShapeInd)%NodeInd(:)
            do j=1,4
              if(this%NodeInd(j)==lPTet(1).or.this%NodeInd(j)==lPTet(2).or.&
              &  this%NodeInd(j)==lPTet(3).or.this%NodeInd(j)==lPTet(4))then
                maskTet(j)=.true.
              end if
            end do
  
negative.code:     &                    Node(this%NodeInd(3))%Pos(:),Node(this%NodeInd(4))%Pos(:))
    end function
    
    !-------------------------------------------------------------
    ! get the k_th neighbour element or facet of this tetrahedron
    !-------------------------------------------------------------
    ! Note: if facet are found, the function will return the negative value of the facet index.
    ! Note: we are getting neighbour element or facet, not neighbour hexahedron or quadrilateral
    !       index, though the neighbour element (facet) is likely to be a tetrahedron (triangle).
    function getTetNeib(this,k)
      class(typeTet),intent(in)::this
      integer,intent(in)::k
      integer getTetNeib,lPTet(4),lPTri(3)
      logical maskTet(4)
      getTetNeib=0
      if(k<1.or.k>4)then
  