voisin31.f

Go to the documentation of this file.

!                   *******************
                    SUBROUTINE voisin31
!                   *******************
!
     &(ifabor,nelem,nelmax,ielm,ikle,sizikl,npoin,nbor,nptfr,
     & lihbor,klog,indpu,iklestr,neleb2)
!
!***********************************************************************
! BIEF   V6P3                                   21/08/2010
!***********************************************************************
!
!brief    BUILDS THE ARRAY IFABOR, WHERE IFABOR(IELEM, IFACE) IS
!+                THE GLOBAL NUMBER OF THE NEIGHBOUR OF SIDE IFACE OF
!+                ELEMENT IELEM (IF THIS NEIGHBOUR EXISTS) AND 0 IF THE
!+                SIDE IS ON THE DOMAIN BOUNDARY.
!
!history  REGINA NEBAUER (LNHE)
!+        22/01/08
!+        V5P8
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (LNHE)
!+        10/09/2012
!+        V6P3
!+   Treatment of prisms cut into tetrahedra added (IELM=51)
!+   and simplifications.
!
!history  S.E.BOURBAN (HRW)
!+        21/03/2017
!+        V7P3
!+   Replacement of the DATA declarations by the PARAMETER associates
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IELM           |-->| 31: TETRAHEDRA
!|                |   | 51: PRISMS CUT INTO TETRAHEDRA
!| IFABOR         |-->| ELEMENTS BEHIND THE EDGES OF A TRIANGLE
!|                |   | IF NEGATIVE OR ZERO, THE EDGE IS A LIQUID
!|                |   | BOUNDARY
!| IKLE           |-->| CONNECTIVITY TABLE.
!| INDPU          |-->| IF NOT 0, INTERFACE POINT IN PARALLEL
!| IKLESTR        |-->| CONNECTIVITY TABLE OF BOUNDARY TRIANGLES
!| KLOG           |-->| CONVENTION FOR SOLID BOUNDARY
!| LIHBOR         |-->| TYPE OF BOUNDARY CONDITIONS ON DEPTH
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELEB2         |-->| NUMBER OF BOUNDARY TRIANGLES (oversized)
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| SIZIKL         |-->| FIRST DIMENSION OF IKLE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_voisin31 => voisin31
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)   :: IELM,NPTFR,NELEM,NELMAX,NPOIN,SIZIKL,KLOG
      INTEGER, INTENT(IN)   :: NBOR(nptfr)
      INTEGER, INTENT(INOUT):: IFABOR(nelmax,4)
      INTEGER, INTENT(IN)   :: IKLE(sizikl,4),LIHBOR(*)
      INTEGER, INTENT(IN)   :: INDPU(*)
      INTEGER, INTENT(IN)   :: NELEB2
      INTEGER, INTENT(IN)   :: IKLESTR(neleb2,3)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
! LOCAL VARIABLES
!
      ! ARRAY WHICH IS THE REVERSE OF NBOR
      ! (GIVES FOR EACH NODE IN THE DOMAIN THE BOUNDARY NODE NUMBER E,
      ! OR 0 IF IT IS AN INTERIOR NODE)
      INTEGER, DIMENSION(:  ), ALLOCATABLE :: NBOR_INV
      ! ARRAY DEFINING THE NUMBER OF ELEMENT (TETRAHEDRONS) NEIGHBOURING
      ! A NODE
      INTEGER, DIMENSION(:), ALLOCATABLE :: NVOIS
      ! ARRAY DEFINING THE IDENTIFIERS OF THE ELEMENTS NEIGHBOURING
      ! EACH NODE
      INTEGER, DIMENSION(:  ), ALLOCATABLE :: NEIGH
      INTEGER, DIMENSION(:,:), ALLOCATABLE :: IKLE_TRI
      INTEGER, DIMENSION(:,:), ALLOCATABLE :: VOIS_TRI
      ! ARRAY DEFINING THE ADDRESSES OF THE VARIOUS ENTRIES IN
      ! ARRAY NEIGH
      INTEGER, DIMENSION(NPOIN)            :: IADR
      ! THE VALUE OF AN ENTRY IN THIS ARRAY
      INTEGER                              :: ADR
!
      ! THE MAXIMUM NUMBER OF ELEMENTS NEIGHBOURING A NODE
      INTEGER :: NMXVOISIN
      INTEGER :: IMAX       ! SIZE OF ARRAY IADR
!
      INTEGER :: NFACE      ! NUMBER OF SIDES TO ELEMENT (TETRA: 4)
      INTEGER :: NBTRI      ! NUMBER OF DEFINED TRIANGLES
!
      INTEGER :: IELEM      ! ELEMENTS COUNTER
      INTEGER :: IELEM2     ! ELEMENTS COUNTER
      INTEGER :: IPOIN      ! DOMAIN NODES COUNTER
      INTEGER :: INOEUD     ! TETRAHEDRONS/TRIANGLES NODES COUNTER
      INTEGER :: IFACE      ! SIDE COUNTER
      INTEGER :: IFACE2     ! SIDE COUNTER
      INTEGER :: ITRI       ! TRIANGLES COUNTER
      INTEGER :: IVOIS      ! NEIGHBOURS COUNTER
      INTEGER :: NV         ! NUMBER OF NEIGHBOURS
!
      INTEGER :: ERR        ! ERROR CODE FOR MEMORY ALLOCATION
!
      LOGICAL :: FOUND      ! FOUND OR NOT ...
!
      INTEGER :: I1, I2, I3 ! THE THREE NODES OF A TRIANGLE
      INTEGER :: M1, M2, M3 ! SAME THING, ORDERED
!
      INTEGER :: I,J,K      ! USEFUL ...
!
      INTEGER :: IR4,IR5,IR6,COMPT
      LOGICAL :: BORD
!
!     DEFINES THE FOUR TRIANGLES OF THE TETRAHEDRON: THE FIRST
!     DIMENSION IS THE NUMBER OF THE TRIANGLE, THE SECOND GIVES
!     THE NODE NUMBERS OF THE NODES OF TETRAHEDRONS WHICH DEFINE IT.
      INTEGER :: SOMFAC(3,4)
      parameter( somfac = reshape( (/
     &    1,2,3 , 4,1,2 , 2,3,4 , 3,4,1 /), shape=(/ 3,4 /) ) )
!
!-----------------------------------------------------------------------
! START
!-----------------------------------------------------------------------
!
! CHECKS FIRST THAT WE ARE INDEED DEALING WITH TETRAHEDRONS. IF NOT,
! GOODBYE.
!
      IF(ielm.EQ.31.OR.ielm.EQ.51) THEN
        nface = 4
      ELSE
        WRITE(lu,99) ielm
99    FORMAT(1x,'VOISIN31: IELM=',1i6,' TYPE OF ELEMENT NOT AVAILABLE')
        CALL plante(1)
        stop
      ENDIF
!
! ALLOCATES TEMPORARY ARRAYS
!
      ALLOCATE(nbor_inv(npoin),stat=err)
      CALL check_allocate(err,'VOISIN31:NBOR_INV')
!
! ALLOCATES TEMPORARY ARRAYS
!
      ALLOCATE(nvois(npoin),stat=err)
      CALL check_allocate(err,'VOISIN31:NVOIS')
!
!-----------------------------------------------------------------------
! STEP 1: COUNTS THE NUMBER OF ELEMENTS NEIGHBOURING A NODE
!-----------------------------------------------------------------------
! COMPUTES THE NUMBER OF ELEMENTS NEIGHBOURING EACH NODE OF THE MESH.
! RESULT: NVOIS(INOEUD) GIVES THE NUMBER OF ELEMENTS NEIGHBOURING
! NODE INOEUD
!
      ! INITIALISES THE NEIGHBOURING ELEMENT COUNTER TO 0
      !
      DO i = 1,npoin
        nvois(i) = 0
      ENDDO
      ! COUNTS THE NEIGHBOURING ELEMENTS
      ! USING THE CONNECTIVITY TABLE, THE COUNTER IS INCREMENTED
      ! EACH TIME THAT AN ELEMENT REFERENCES NODE IPOIN
      ! LOOP ON THE 4 NODES OF THE ELEMENT
      DO inoeud = 1, 4
        ! LOOP ON THE ELEMENTS
        DO ielem = 1,nelem
          ! THE ID OF NODE I OF ELEMENT IELEM
          ipoin        = ikle( ielem , inoeud )
          ! INCREMENT THE COUNTER
          nvois(ipoin) = nvois(ipoin) + 1
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
! STEP 2: DETERMINES THE SIZE OF ARRAY NEIGH() AND OF AUXILIARY
! ARRAY TO INDEX NEIGH. ALLOCATES NEIGH
!-----------------------------------------------------------------------
! CREATES AN ARRAY WHICH WILL CONTAIN THE IDENTIFIERS OF THE ELEMENTS
! NEIGHBOURING EACH NODE. SINCE THE NUMBER OF NEIGHBOURS IS A PRIORI
! DIFFERENT FOR EACH NODE, AND IN AN EFFORT NOT TO CREATE A (TOO) BIG
! ARRAY FOR NO REASON, AN AUXILIARY ARRAY IS REQUIRED THAT GIVES THE
! ADDRESS OF THE ENTRIES FOR A GIVEN NODE. THIS ARRAY HAS AS MANY
! ENTRIES AS THERE ARE NODES.
! WILL ALSO COMPUTE THE MAXIMUM NUMBER OF NEIGHBOURS, SOME TIME.
!
      ! THE FIRST ENTRY IN THE ID OF THE NEIGHBOURS ARRAY IS 1
      adr       = 1
      iadr(1)   = adr
      ! THE MAX NUMBER OF NEIGHBOURING ELEMENTS
      nv        = nvois(1)
      nmxvoisin = nv
!
      DO ipoin = 2,npoin
        ! ADDRESS FOR THE OTHER ENTRIES:
        adr         = adr + nv
        iadr(ipoin) = adr
        nv          = nvois(ipoin)
        ! IDENTIFIES THE MAX. NUMBER OF NEIGHBOURS
        nmxvoisin   = max(nmxvoisin,nv)
      ENDDO
!
      ! THE TOTAL NUMBER OF NEIGHBOURING ELEMENTS FOR ALL THE NODES
      ! GIVES THE SIZE OF THE NEIGHBOURS ARRAY:
      imax = iadr(npoin) + nvois(npoin)
      ! ALLOCATES THE ARRAY CONTAINING THE IDENTIFIERS OF THE ELEMENTS
      ! NEIGHBOURING EACH NODE
      ALLOCATE(neigh(imax),stat=err)
      IF(err.NE.0) GOTO 999
!
!-----------------------------------------------------------------------
! STEP 3: INITIALISES NEIGH
!-----------------------------------------------------------------------
! NEEDS TO FILL NEIGH NOW THAT IT'S BEEN ALLOCATED: I.E.
! STARTS AGAIN THE LOOP ON THE 4 NODES OF EACH ELEMENT AND THIS
! TIME, ALSO STORES THE IDENTIFIER IN ARRAY NEIGH.
!
!
      ! RE-INITIALISES THE COUNTER OF THE NEIGHBOURING ELEMENTS TO 0,
      ! TO KNOW WHERE WE ARE AT
      DO i = 1,npoin
        nvois(i) = 0
      ENDDO
!     NVOIS(:) = 0
      ! FOR EACH NODE OF THE ELEMENTS, STORES THE IDENTIFIER
      DO inoeud = 1, 4  ! LOOP ON THE ELEMENT NODES
        DO ielem=1,nelem ! LOOP ON THE ELEMENTS
          ipoin     = ikle( ielem , inoeud )
          ! ONE MORE NEIGHBOUR
          nv           = nvois(ipoin) + 1
          nvois(ipoin) = nv
          ! STORES THE IDENTIFIER OF THE NEIGHBOURING ELEMENT IN THE ARRAY
          neigh(iadr(ipoin)+nv) = ielem
        ENDDO ! END OF LOOP ELEMENTS
      ENDDO  ! END OF LOOP NODES
!
!-----------------------------------------------------------------------
! STEP 4: IDENTIFIES COMMON FACES OF THE TETRAHEDRONS AND FILLS IN
! ARRAY IFABOR
!-----------------------------------------------------------------------
! TO IDENTIFY FACES COMMON TO THE ELEMENTS :
! FROM THE ELEMENTS THAT SHARE A NODE, AT LEAST 2 SHARE A FACE
! (IF THE NODE IS NOT A BOUNDARY NODE).
! THE ALGORITHM PRINCIPLE:
! BUILDS THE TRIANGLES OF THE TETRAHEDRON FACES, ONCE HAVE IDENTIFIED
! THOSE THAT SHARE NODE IPOIN.
! IF 2 TRIANGLES SHARE THE SAME NODES, IT MEANS THAT THE TETRAHEDRONS
! DEFINING THEM ARE NEIGHBOURS.
! IF NO NEIGHBOUR CAN BE FOUND, IT MEANS THAT THE TRIANGLE IS A
! BOUNDARY FACE.
! BASED ON THE ASSUMPTION THAT A TRIANGLE CANNOT BE DEFINED BY MORE
! THAN 2 TETRAHEDRONS.
! IF THAT WAS NOT THE CASE, IT WOULD MEAN THAT THERE WAS A PROBLEM WITH
! THE MESH; AND THIS IS NOT CATERED FOR ...
!
! ADVANTAGES:
! SAVES QUITE A BIT OF MEMORY, BY STORING THE TRIANGLES AROUND A NODE.
! DISADVANTAGES:
! COULD BE DOING TOO MANY (USELESS) COMPUTATIONS (TO GET TO THE STAGE
! WHERE THE CONNECTIVITY TABLE FOR THE TRIANGLES IS DEFINED)
! COULD MAYBE SKIP THIS STEP BY CHECKING IF IFABOR ALREADY CONTAINS
! SOMETHING OR NOT ...
!
! BUILDS THE CONNECTIVITY TABLE FOR THE TRIANGLES
! THIS CONNECTIVITY TABLE IS NOT SUPPOSED TO MAKE A LIST OF ALL
! THE TRIANGLES, BUT MERELY THOSE AROUND A NODE.
! THE MAXIMUM NUMBER OF (TETRAHEDRONS) NEIGHBOURS IS KNOWN FOR A
! NODE. IN THE WORST CASE, THE NODE IS A BOUNDARY NODE.
! WILL MAXIMISE (A LOT) BY ASSUMING THAT THE MAXIMUM NUMBER OF
! TRIANGLES AROUND A NODE CAN BE THE NUMBER OF NEIGHBOURING
! TETRAHEDRONS.
! ALSO BUILDS THE ARRAY VOIS_TRI CONTAINING THE ID OF THE TETRAHEDRON
! ELEMENT THAT DEFINED IT FIRST (AND THAT WILL BE NEIGHBOUR TO THE
! TETRAHEDRON THAT WILL FIND THAT ANOTHER ONE ALREADY DEFINES IT)
! THIS ARRAY HAS 2 ENTRIES : THE ID OF THE ELEMENT AND THE ID OF THE SIDE.
!
      nbtri = nmxvoisin * 3
!
      ALLOCATE(ikle_tri(nbtri,3),stat=err)
      IF(err.NE.0) GOTO 999
      ALLOCATE(vois_tri(nbtri,2),stat=err)
      IF(err.NE.0) GOTO 999
!
!     ALL FACES PRESUMED SOLID TO START WITH...
!
      ! ASSUMING SOLID BOUNDARIES EVERYWHERE
      ifabor(:,:) = -1
!
      ! LOOP ON ALL THE NODES IN THE MESH
      DO ipoin = 1, npoin
        ! FOR EACH NODE, CHECKS THE NEIGHBOURING TETRAHEDRON ELEMENTS
        ! (MORE PRECISELY: THE TRIANGULAR FACES THAT MAKE IT)
        ! RE-INITIALISES THE CONNECTIVITY TABLE FOR THE TETRAHEDRON
        ! TRIANGLES TO 0, AND THE NUMBER OF TRIANGLES WHICH HAVE BEEN
        ! FOUND:
        !
        ikle_tri(:,:) = 0
        ! SAME THING FOR THE ARRAY THAT IDENTIFIES WHICH ELEMENT HAS
        ! ALREADY DEFINED THE TRIANGLE :
        vois_tri(:,:) = 0
        ! STARTS COUNTING THE TRIANGLES AGAIN:
        nbtri         = 0
        nv            = nvois(ipoin)
        adr           = iadr(ipoin)
        DO ivois = 1, nv
          ! THE IDENTIFIER OF THE NEIGHBOURING ELEMENT IVOIS TO
          ! THE NODE IPOIN:
          ielem = neigh(adr+ivois)
          ! LOOP ON THE 4 SIDES OF THIS ELEMENT
          DO iface = 1 , nface
            ! IF THIS SIDE ALREADY HAS A NEIGHBOUR, THAT'S
            ! ENOUGH AND GOES TO NEXT.
            ! OTHERWISE, LOOKS FOR IT...
            IF ( ifabor(ielem,iface) .LE. 0 ) THEN
            ! EACH FACE DEFINES A TRIANGLE. THE TRIANGLE IS
            ! GIVEN BY 3 NODES.
            i1 = ikle(ielem,somfac(1,iface))
            i2 = ikle(ielem,somfac(2,iface))
            i3 = ikle(ielem,somfac(3,iface))
            ! THESE 3 NODES ARE ORDERED, M1 IS THE NODE WITH
            ! THE SMALLEST IDENTIFIER, M3 THAT WITH THE
            ! LARGEST IDENTIFIER AND M2 IS IN THE MIDDLE:
            m1 = max(i1,(max(i2,i3)))
            m3 = min(i1,(min(i2,i3)))
            m2 = i1+i2+i3-m1-m3
            ! GOES THROUGH THE ARRAY WITH TRIANGLES ALREADY DEFINED
            ! TO SEE IF ONE OF THEM BEGINS WITH M1.
            ! IF THAT'S THE CASE, CHECKS THAT IT ALSO HAS NODES
            ! M2 AND M3. IF THAT'S THE CASE, HAS FOUND A NEIGHBOUR.
            ! OTHERWISE, A NEW TRIANGLE ENTRY IS CREATED.
            !
            found = .false.
            DO itri = 1, nbtri
              IF ( ikle_tri(itri,1) .EQ. m1 ) THEN
                IF ( ikle_tri(itri,2) .EQ. m2 .AND.
     &            ikle_tri(itri,3) .EQ. m3 ) THEN
                  ! FOUND IT! ALL IS WELL.
                  ! STORES THE INFORMATION IN VOIS_TRI.
                  ! (I.E. THE ELEMENT THAT HAS ALREADY
                  ! DEFINED THE TRIANGLE AND THE FACE)
                  ielem2 = vois_tri(itri,1)
                  iface2 = vois_tri(itri,2)
                  IF ( ielem2 .EQ. ielem ) THEN
                    WRITE(lu,909) ielm
909                 FORMAT(1x,'VOISIN: IELM=',1i6,',
     &              NEIGHBOUR PROBLEM')
                    CALL plante(1)
                    stop
                  END IF
                  ! TO BE SURE :
                  IF ( ielem2 .EQ. 0 .OR.
     &                 iface2 .EQ. 0 ) THEN
                    WRITE(lu,919) ielem2,iface2
919                 FORMAT(1x,'VOISIN31:UNDEFINED TRIANGLE,
     &              IELEM=',1i6,'IFACE=',1i6)
                    CALL plante(1)
                    stop
                  END IF
                  ! THE ELEMENT AND ITS NEIGHBOUR : STORES
                  ! THE CONNECTION IN IFABOR.
                  ifabor(ielem ,iface ) = ielem2
                  ifabor(ielem2,iface2) = ielem
                  found = .true.
                END IF
              END IF
            END DO
            ! NO, THIS TRIANGLE WAS NOT ALREADY THERE; THEREFORE
            ! CREATES A NEW ENTRY.
            IF ( .NOT. found) THEN
                nbtri             = nbtri + 1
                ikle_tri(nbtri,1) = m1
                ikle_tri(nbtri,2) = m2
                ikle_tri(nbtri,3) = m3
                vois_tri(nbtri,1) = ielem
                vois_tri(nbtri,2) = iface
            END IF
          END IF ! IFABOR 0
          END DO ! END OF LOOP ON FACES OF THE NEIGHBOURING ELEMENTS
!
        END DO ! END OF LOOP ON ELEMENTS NEIGHBOURING THE NODE
      END DO ! END OF LOOP ON NODES
!
      DEALLOCATE(neigh)
      DEALLOCATE(ikle_tri)
      DEALLOCATE(vois_tri)
!
!-----------------------------------------------------------------------
!     STEP 5: FACES BETWEEN DIFFERENT COMPUTATION DOMAIN
!-----------------------------------------------------------------------
!
      IF(ncsize.GT.1) THEN
        !
        ! WHY A DIFFERENT ALGORITHM ?
        ! a) Partitionning method : 51 submeshes are obtained from a 2D mesh
        ! b) No iklestr is available
        ! c) Some pathological cases with 31
        IF (ielm.EQ.51) THEN
          !
          DO iface=1,nface
            DO ielem=1,nelem
!     IF A FACE HAS 3 POINTS WHICH ARE INTERFACES BETWEEN SUB-DOMAINS
!     IT IS ASSIGNED A VALUE OF -2
!
              i1=ikle(ielem,somfac(1,iface))
              i2=ikle(ielem,somfac(2,iface))
              i3=ikle(ielem,somfac(3,iface))
!
              IF( indpu(i1).NE.0.AND.
     &            indpu(i2).NE.0.AND.
     &            indpu(i3).NE.0     ) ifabor(ielem,iface)=-2
!
            ENDDO
          ENDDO
!
        ELSE IF (ielm.EQ.31) THEN
!
          DO iface=1,nface
            DO ielem=1,nelem
!     IF A FACE HAS 3 POINTS WHICH ARE INTERFACES BETWEEN SUB-DOMAINS
!     IT IS ASSIGNED A VALUE OF -2

!     PATHOLOGIC CASES (THANKS TO OLIVIER BOITEAU, EDF R&D/SINETICS)
!     a) 3 nodes may have INDPU /= 0, this doesn't mean automatically
!     that the triangle is an interface one (case of true unstructured tetra mesh)
!     b) a border triangle may have INDPU /= 0 on its 3 nodes
!     => priority to border triangles
!     c) 3 nodes of a triangle may be on the boundary and have INDPU /= 0
!     this doesn't mean automatically that the triangle is on the boundary.

              IF (ifabor(ielem,iface).EQ.-1) THEN
!
                i1=ikle(ielem,somfac(1,iface))
                i2=ikle(ielem,somfac(2,iface))
                i3=ikle(ielem,somfac(3,iface))

                IF ( indpu(i1).NE.0.AND.
     &               indpu(i2).NE.0.AND.
     &               indpu(i3).NE.0     ) THEN

!     THESE ARE INTERFACE NODES; DO THEY CORRESPOND TO A
!     (VIRTUAL) INTERFACE TRIANGLE OR TO A BOUNDARY TRIANGLE ?

!     THE FOLLOWING TEST IS NOT SUFFICIENT
!     IF (.NOT. (IR5.EQ.1.AND.IR4.EQ.1.AND.IR6.EQ.1)) IFABOR(IELEM,IFACE)=-2
!
                  bord=.false.
                  ir4=0
                  ir5=0
                  ir6=0
                  DO j=1,nptfr
                    IF (i1.EQ.nbor(j)) ir5=1
                    IF (i2.EQ.nbor(j)) ir4=1
                    IF (i3.EQ.nbor(j)) ir6=1
                  ENDDO ! J

!     THEY ARE ALSO BOUNDARY NODES
                  IF (ir5.EQ.1.AND.ir4.EQ.1.AND.ir6.EQ.1) THEN
!
                    DO j=1,neleb2
!     IT IS A BOUNDARY TRIANGLE
                      compt=0
                      DO i=1,3
                        IF (iklestr(j,i)==i1) compt=compt+1
                        IF (iklestr(j,i)==i2) compt=compt+10
                        IF (iklestr(j,i)==i3) compt=compt+100
                      ENDDO
!     THESE 3 NODES INDEED BELONG TO THE SAME BOUNDARY TRIANGLE
                      IF (compt==111) THEN
                        bord=.true.
!     WRITE(LU,*)'VERTICES OF A BOUNDARY TRIANGLE'
                        EXIT
                      ENDIF
                    ENDDO
                  ENDIF
                  IF (.NOT.bord) THEN
!     THESE 3 NODES BELONG TO AN INTERFACE MESH
!     WRITE(LU,*) 'INTERFACE NODES'
                    ifabor(ielem,iface)=-2
                  ENDIF
                ENDIF
              ENDIF
            ENDDO
          ENDDO
        ENDIF
      ENDIF
!
!
!-----------------------------------------------------------------------
!
!     WITH ELEMENT 51 THE DIFFERENCE BETWEEN SOLID AND LIQUID
!     BOUNDARIES IS DONE LATER
!     AND I DOUBT THAT THE CODE BELOW WITH LIHBOR CAN WORK WITH
!     REAL TETRAHEDRA
!
      IF((ielm.EQ.51).OR.(ielm.EQ.31)) GO TO 1000
!
!-----------------------------------------------------------------------
!
!  IFABOR DISTINGUISHES BETWEEN THE BOUNDARY FACES AND THE LIQUID FACES
!
!  INITIALISES NBOR_INV TO 0
!
      DO ipoin=1,npoin
        nbor_inv(ipoin) = 0
      ENDDO
!
!  CONNECTS GLOBAL NUMBERING TO BOUNDARY NUMBERING
!
      DO k = 1, nptfr
        nbor_inv(nbor(k)) = k
      ENDDO
!
!  LOOP ON ALL THE SIDES OF ALL THE ELEMENTS:
!
      DO iface = 1 , nface
      DO ielem = 1 , nelem
!
      IF(ifabor(ielem,iface).EQ.-1) THEN
!
!       IT IS A TRUE BOUNDARY SIDE (IN PARALLEL MODE THE INTERNAL SIDES
!                                   ARE INDICATED WITH -2).
!       GLOBAL NUMBERS OF THE NODES OF THE SIDE :
!
        i1 = ikle( ielem , somfac(1,iface) )
        i2 = ikle( ielem , somfac(2,iface) )
        i3 = ikle( ielem , somfac(3,iface) )
!
!       A LIQUID SIDE IS IDENTIFIED WITH THE BOUNDARY CONDITION ON H
!
        IF(nptfr.GT.0) THEN
        IF(lihbor(nbor_inv(i1)).NE.klog.AND.lihbor(nbor_inv(i2)).NE.klog
     &      .AND.lihbor(nbor_inv(i3)).NE.klog  ) THEN
!         LIQUID SIDE : IFABOR=0  SOLID SIDE : IFABOR=-1
          ifabor(ielem,iface)=0
        ENDIF
        ENDIF
!
      ENDIF
!
      ENDDO ! IELEM
      ENDDO ! IFACE
!
!
1000  CONTINUE
!
!-----------------------------------------------------------------------
!
      RETURN
!
!-----------------------------------------------------------------------
!
999   WRITE(lu,2000) err
2000  FORMAT(1x,'VOISIN31: ERROR DURING ALLOCATION OF MEMORY: ',/,1x,
     &            'ERROR CODE: ',1i6)
      CALL plante(1)
      stop
!
!-----------------------------------------------------------------------
!
      END