voisin.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE voisin
!                   *****************
!
     &(ifabor,nelem,nelmax,ielm,ikle,sizikl,
     & npoin,nachb,nbor,nptfr,iadr,nvois)
!
!***********************************************************************
! BIEF   V6P1                                   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  OLIVIER BOITEAU (SINETICS)
!+        19/02/08
!+
!+   SIZE OF NACHB
!
!history  J-M HERVOUET (LNHE)
!+        16/06/08
!+        V5P9
!+   MODIFICATION TO THE MAX NUMBER OF NEIGHBOURS
!
!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  S.E.BOURBAN (HRW)
!+        21/03/2017
!+        V7P3
!+   Replacement of the DATA declarations by the PARAMETER associates
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IADR           |<->| WORK ARRAY
!| IELM           |-->| 11: TRIANGLES
!|                |   | 21: QUADRILATERALS
!| IFABOR         |-->| ELEMENTS BEHIND THE EDGES OF A TRIANGLE
!|                |   | IF NEGATIVE OR ZERO, THE EDGE IS A LIQUID
!|                |   | BOUNDARY
!| IKLE           |-->| CONNECTIVITY TABLE.
!| NACHB          |-->| SUB-DOMAINS NEIGHBOURS OF INTERFACE POINTS IN PARALLEL
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NVOIS          |<--| NUMBER OF NEIGHBOURS OF POINTS
!| SIZIKL         |-->| FIRST DIMENSION OF IKLE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_voisin => voisin
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)    :: NPTFR,SIZIKL,NELEM,NELMAX,IELM,NPOIN
      INTEGER, INTENT(IN)    :: NBOR(nptfr),NACHB(nbmaxnshare,nptir)
      INTEGER, INTENT(IN)    :: IKLE(sizikl,*)
      INTEGER, INTENT(INOUT) :: IFABOR(nelmax,*)
      INTEGER, INTENT(INOUT) :: NVOIS(npoin),IADR(npoin)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER NFACE,NDP,KEL,IMAX,IFACE,IELEM,M1,M2,IV,IELEM2,IFACE2
      INTEGER I,J,ERR,IR1,IR2,IR3,IR4,I1,I2,IDIMAT
!
      INTEGER :: SOMFAC(2,4,2)
      parameter( somfac = reshape( (/
     &      1,2 , 2,3 , 3,1 , 0,0   ,
     &      1,2 , 2,3 , 3,4 , 4,1 /), shape=(/ 2,4,2 /) ) )
!
!     DYNAMICALLY ALLOCATES THE WORKING ARRAYS
!
      INTEGER, DIMENSION(:), ALLOCATABLE :: MAT1,MAT2,MAT3
!
!-----------------------------------------------------------------------
!
      IF(ielm.EQ.21) THEN
!       QUADRILATERALS
        nface = 4
!       NUMBER OF POINTS PER ELEMENT
        ndp = 4
!       ADDRESS IN SOMFAC
        kel = 2
      ELSEIF(ielm.EQ.11.OR.ielm.EQ.41.OR.ielm.EQ.51) THEN
!       TRIANGLES
        nface = 3
!       NUMBER OF POINTS PER ELEMENT
        ndp = 3
!       ADDRESS IN SOMFAC
        kel = 1
      ELSE
        WRITE(lu,99) ielm
99      FORMAT(1x,'VOISIN: IELM=',1i6,' TYPE OF ELEMENT NOT AVAILABLE')
        CALL plante(1)
        stop
      ENDIF
!
!     IDIMAT IS BIGGER THAN THE SUM OF THE NUMBER OF NEIGHBOURS OF
!     ALL THE POINTS (NEIGHBOUR = CONNECTED BY A SEGMENT)
!
      idimat = ndp*2*nelem
!
      ALLOCATE(mat1(idimat),stat=err)
      ALLOCATE(mat2(idimat),stat=err)
      ALLOCATE(mat3(idimat),stat=err)
!
      IF(err.NE.0) THEN
        WRITE(lu,2000) err
2000    FORMAT(1x,'VOISIN: ERROR DURING ALLOCATION OF MEMORY: ',/,1x,
     &            'ERROR CODE: ',1i6)
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!  ARRAY NVOIS FOR EACH POINT
!  BEWARE : NVOIS IS BIGGER THAN THE ACTUAL NUMBER OF NEIGHBOURS
!           THE SUM OF NVOIS WILL GIVE IDIMAT
!
      DO i=1,npoin
        nvois(i) = 0
      ENDDO
!
      DO iface = 1,nface
        DO ielem=1,nelem
          i1 = ikle( ielem , somfac(1,iface,kel) )
          i2 = ikle( ielem , somfac(2,iface,kel) )
          nvois(i1) = nvois(i1) + 1
          nvois(i2) = nvois(i2) + 1
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
!  ADDRESSES OF EACH POINT IN A STRUCTURE OF TYPE COMPACT MATRIX
!
!
      iadr(1) = 1
      DO i= 2,npoin
        iadr(i) = iadr(i-1) + nvois(i-1)
      ENDDO ! I
!
      imax = iadr(npoin) + nvois(npoin) - 1
      IF(imax.GT.idimat) THEN
        WRITE(lu,52) idimat,imax
52      FORMAT(1x,'VOISIN: SIZE OF MAT1,2,3 (',1i9,') TOO SHORT',/,
     &         1x,'MINIMUM SIZE: ',1i9)
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!  INITIALISES THE COMPACT MATRIX TO 0
!
      DO i=1,imax
        mat1(i) = 0
      ENDDO
!
!-----------------------------------------------------------------------
!
!  LOOP ON THE SIDES OF EACH ELEMENT:
!
      DO iface = 1 , nface
      DO ielem = 1 , nelem
!
      ifabor(ielem,iface) = -1
!
!       GLOBAL NODE NUMBERS FOR THE SIDE:
!
        i1 = ikle( ielem , somfac(1,iface,kel) )
        i2 = ikle( ielem , somfac(2,iface,kel) )
!
!       ORDERED GLOBAL NUMBERS:
!
        m1 = min0(i1,i2)
        m2 = max0(i1,i2)
!
        DO iv = 1,nvois(m1)
!
          IF(mat1(iadr(m1)+iv-1).EQ.0) THEN
            mat1(iadr(m1)+iv-1)=m2
            mat2(iadr(m1)+iv-1)=ielem
            mat3(iadr(m1)+iv-1)=iface
            GO TO 81
          ELSEIF(mat1(iadr(m1)+iv-1).EQ.m2) THEN
            ielem2 = mat2(iadr(m1)+iv-1)
            iface2 = mat3(iadr(m1)+iv-1)
            ifabor(ielem,iface) = ielem2
            ifabor(ielem2,iface2) = ielem
            GO TO 81
          ENDIF
!
        ENDDO ! IV
!
        WRITE(lu,83)
83      FORMAT(1x,'VOISIN : ERROR IN THE MESH             ',/,1x,
     &            '         MAYBE SUPERIMPOSED POINTS     ')
        CALL plante(1)
        stop
!
81      CONTINUE
!
      ENDDO ! IELEM
      ENDDO ! IFACE
!
!  COULD TRY SOMETHING A BIT LIGHTER
!  USING INDPU FOR EXAMPLE
!
      IF(ncsize.GT.1) THEN
!
      DO iface=1,nface
      DO ielem=1,nelem
!
!  SOME BOUNDARY SIDES ARE INTERFACES BETWEEN SUB-DOMAINS IN
!  ACTUAL FACT: THEY ARE ASSIGNED A VALUE -2 INSTEAD OF -1
!
      IF(ifabor(ielem,iface).EQ.-1) THEN
!
        i1 = ikle( ielem , somfac(1,iface,kel) )
        i2 = ikle( ielem , somfac(2,iface,kel) )
!
        ir1=0
        ir2=0
!
        IF(nptir.GT.0) THEN
          DO j=1,nptir
            IF(i1.EQ.nachb(1,j)) ir1=1
            IF(i2.EQ.nachb(1,j)) ir2=1
          ENDDO ! J
        ENDIF
!
        IF(ir1.EQ.1.AND.ir2.EQ.1) THEN
!         INTERFACE SEGMENT DETECTED, CHECKS WHETHER IT IS NOT
!         ALSO A TRUE BOUNDARY SIDE
          ir3=0
          ir4=0
          DO j=1,nptfr
            IF(i1.EQ.nbor(j)) ir3=1
            IF(i2.EQ.nbor(j)) ir4=1
          ENDDO ! J
!         PRIORITY TO THE TRUE BOUNDARY SIDES
          IF(ir3.EQ.0.OR.ir4.EQ.0) THEN
            ifabor(ielem,iface)=-2
          ENDIF
        ENDIF
!
      ENDIF
!
      ENDDO ! IELEM
      ENDDO ! IFACE
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      DEALLOCATE(mat1)
      DEALLOCATE(mat2)
      DEALLOCATE(mat3)
!
!-----------------------------------------------------------------------
!
      RETURN
      END