elebd.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE elebd
!                   ****************
!
     &(nelbor,nulone,kp1bor,ifabor,nbor,ikle,sizikl,iklbor,nelem,nelmax,
     & npoin,nptfr,ielm,lihbor,klog,iseg,t1,t2,t3,nelebx,
     & neleb)
!
!***********************************************************************
! BIEF   V7P2
!***********************************************************************
!
!brief    BUILDING DATA STRUCTURES TO NAVIGATE IN A 2D MESH.
!+
!+            1) ARRAYS NELBOR AND NULONE,
!+
!+            2) ARRAY KP1BOR,
!+
!+            3) DISTINGUISHES IN THE ARRAY IFABOR FOR
!+                   SOLID BOUNDARY FACES OR LIQUID FACES,
!+
!+            4) IKLBOR, CONNECTIVITY TABLE OF BOUNDARY ELEMENTS
!
!history  J-M HERVOUET (LNHE)
!+        23/06/2008
!+        V5P9
!+
!
!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)
!+        18/03/2013
!+        V6P3
!+   Error messages more accurate.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        13/03/2014
!+        V7P0
!+   Now written to enable different numbering of boundary points and
!+   boundary segments.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        16/06/2014
!+        V7P0
!+   New possible errors in the mesh now stopped.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        25/04/2016
!+        V7P2
!+   Checking wrong meshes where a point appears on several different
!+   boundaries.
!
!history  S.E.BOURBAN (HRW)
!+        21/03/2017
!+        V7P3
!+   Replacement of the DATA declarations by the PARAMETER associates
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IELM           |-->| TYPE D'ELEMENT.
!|                |   | 11 : TRIANGLES.
!|                |   | 21 : QUADRILATERES.
!| IFABOR         |-->| TABLEAU DES VOISINS DES FACES.
!| IKLBOR         |<--| CONNECTIVITY TABLE OF BOUNDARY ELEMENTS
!| IKLE           |-->| NUMEROS GLOBAUX DES POINTS DE CHAQUE ELEMENT.
!| ISEG           |---|
!| KLOG           |-->| CONVENTION POUR LA CONDITION LIMITE DE PAROI
!| KP1BOR         |<--| NUMERO DU POINT SUIVANT LE POINT DE BORD K.
!| LIHBOR         |-->| TYPES DE CONDITIONS AUX LIMITES SUR H
!| NBOR           |-->| NUMERO GLOBAL DU POINT DE BORD K.
!| NELBOR         |<--| NUMERO DE L'ELEMENT ADJACENT AU KIEME SEGMENT
!| NELEM          |-->| NOMBRE TOTAL D'ELEMENTS DANS LE MAILLAGE.
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NOMBRE TOTAL DE POINTS DU DOMAINE.
!| NPTFR          |-->| NOMBRE DE POINTS FRONTIERES.
!| NULONE         |<--| NUMERO LOCAL D'UN POINT DE BORD DANS
!|                |   | L'ELEMENT ADJACENT DONNE PAR NELBOR
!| SIZIKL         |-->| FIRST DIMENSION OF IKLE
!| T2             |---|
!| T3             |---|
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_elebd => elebd
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)    :: KLOG,NELMAX,NELEM,SIZIKL,NELEBX
      INTEGER, INTENT(INOUT) :: NELEB
      INTEGER, INTENT(IN)    :: NPOIN,NPTFR,IELM
      INTEGER, INTENT(INOUT) :: NELBOR(nelebx),NULONE(nelebx,2)
      INTEGER, INTENT(INOUT) :: KP1BOR(nptfr,2)
      INTEGER, INTENT(INOUT) :: NBOR(*)
      INTEGER, INTENT(INOUT) :: IFABOR(nelmax,*)
      INTEGER, INTENT(IN)    :: IKLE(sizikl,*)
      INTEGER, INTENT(IN)    :: LIHBOR(nptfr)
      INTEGER, INTENT(INOUT) :: IKLBOR(nelebx,2)
      INTEGER, INTENT(IN)    :: ISEG(*)
      INTEGER, INTENT(INOUT) :: T1(npoin),T2(npoin),T3(npoin)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,NFACE,NPT,KEL,IPOIN
      INTEGER K,IFACE,I1,I2,N1,N2,IPT,IEL,I,K1,K2
!
      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 /) ) )
!
!-----------------------------------------------------------------------
!
      IF(ielm.EQ.11.OR.ielm.EQ.41.OR.ielm.EQ.51) THEN
!       TRIANGLES
        nface = 3
        npt = 3
        kel = 1
      ELSE
        WRITE(lu,901) ielm
901     FORMAT(1x,'ELEBD: IELM=',1i6,' UNKNOWN TYPE OF ELEMENT')
        CALL plante(1)
        stop
      ENDIF
!
!  INITIALISES T1,2,3 TO 0
!
      DO ipoin=1,npoin
        t1(ipoin) = 0
        t2(ipoin) = 0
        t3(ipoin) = 0
      ENDDO
!
!  STORES K IN TRAV(*,3), ADDRESS NBOR(K)
!  GIVES CORRESPONDENCE GLOBAL --> BOUNDARY NUMBER
!  FINDS POINTS WITH MORE THAN ONE BOUNDARY
!
      DO k = 1, nptfr
        t3(nbor(k)) = k
        t1(nbor(k)) = t1(nbor(k)) +1
      ENDDO
      DO i = 1, npoin
        IF(t1(i).GT.1) THEN
          WRITE(lu,*) 'ELEBD: POINT ',i
          WRITE(lu,*) 'APPEARS ON ',t1(i)
          WRITE(lu,*) 'DIFFERENT BOUNDARIES'
          WRITE(lu,*) 'WRONG MESH'
          CALL plante(1)
          stop
        ENDIF
      ENDDO
!
!  REINITIALISES T1 TO 0 (IS THIS USEFUL ??)
!
      DO ipoin=1,npoin
        t1(ipoin) = 0
      ENDDO
!
!  LOOP ON ALL THE FACES OF ALL THE ELEMENTS:
!
      DO iface = 1 , nface
        DO ielem = 1 , nelem
          IF(ifabor(ielem,iface).EQ.-1) THEN
!           THIS IS A TRUE BOUNDARY FACE
!           INTERNAL FACES ARE MARKED WITH -2 IN PARALLELE MODE
!           GLOBAL NUMBERS OF THE FACE POINTS :
            i1 = ikle( ielem , somfac(1,iface,kel) )
            i2 = ikle( ielem , somfac(2,iface,kel) )
!           STORES IN T1 AND T2 (ADDRESS I1) : I2 AND IELEM
            t1(i1) = i2
            t2(i1) = ielem
!           A LIQUID FACE IS RECOGNIZED BY THE BOUNDARY CONDITION ON H
            IF(nptfr.GT.0) THEN
              IF(t3(i1).NE.0.AND.t3(i2).NE.0) THEN
                IF(lihbor(t3(i1)).NE.klog.AND.
     &             lihbor(t3(i2)).NE.klog) THEN
!                 LIQUID FACE : IFABOR=0  SOLID FACE : IFABOR=-1
                  ifabor(ielem,iface)=0
                ENDIF
              ELSE
                IF(t3(i1).EQ.0) THEN
                  WRITE(lu,*) 'POINT ',i1,' IS ON A BOUNDARY'
                ENDIF
                IF(t3(i2).EQ.0) THEN
                  WRITE(lu,*) 'POINT ',i2,' IS ON A BOUNDARY'
                ENDIF
                WRITE(lu,*) 'BUT NOT IN THE LIST OF BOUNDARY POINTS'
                WRITE(lu,*) 'OR THE NUMBER OF LINES IN THE BOUNDARY'
                WRITE(lu,*) 'CONDITIONS FILE IS GREATER THAN:',nptfr
                WRITE(lu,*) 'NUMBER TAKEN IN THE GEOMETRY FILE'
                CALL plante(1)
                stop
              ENDIF
            ENDIF
!
          ENDIF
        ENDDO ! IELEM
      ENDDO ! IFACE
!
!     LOOP ON ALL THE BOUNDARY POINTS TO BUILD KP1BOR:
!
      IF(nptfr.GT.0) THEN
        DO i = 1 , npoin
          IF(t1(i).NE.0) THEN
!           FOLLOWING POINT
            kp1bor(t3(i),1)=t3(t1(i))
!           PRECEDING POINT
            kp1bor(t3(t1(i)),2)=t3(i)
            nelbor(t3(i))=t2(i)
          ENDIF
        ENDDO
      ENDIF
!
!     DUMMY VALUES IN KP1BOR WHEN THE FOLLOWING POINT IS IN ANOTHER SUB-DOMAIN
!     NELBOR AND NULONE SET TO 0
!
      IF(ncsize.GT.1) THEN
        DO k1=1,nptfr
          IF(iseg(k1).GT.0) THEN
            kp1bor(k1,1)=k1
            nelbor(k1)=0
            nulone(k1,1)=0
            nulone(k1,2)=0
          ELSEIF(iseg(k1).EQ.-9999) THEN
            kp1bor(k1,1)=k1
            kp1bor(k1,2)=k1
            nelbor(k1)  =0
            nulone(k1,1)=0
            nulone(k1,2)=0
          ELSEIF(iseg(k1).LT.0) THEN
            kp1bor(k1,2)=k1
          ENDIF
        ENDDO
      ENDIF
!
! COMPUTES ARRAY NULONE
!
      DO k1=1,nptfr
!
      IF(ncsize.GT.1) THEN
        IF(iseg(k1).GT.0.OR.iseg(k1).EQ.-9999) cycle
      ENDIF
!
      k2=kp1bor(k1,1)
      iel = nelbor(k1)
      n1  = nbor(k1)
      IF(k2.LT.0.OR.k2.GT.nptfr) print*,'K2=',k2,' K1=',k1
      n2  = nbor(k2)
!
      i1 = 0
      i2 = 0
!
      DO ipt=1,npt
        IF(ikle(iel,ipt).EQ.n1) THEN
          nulone(k1,1) = ipt
          i1 = 1
        ENDIF
        IF(ikle(iel,ipt).EQ.n2) THEN
          nulone(k1,2) = ipt
          i2 = 1
        ENDIF
      ENDDO
!
      IF(i1.EQ.0.OR.i2.EQ.0) THEN
        WRITE(lu,811) ipt,k1,iel
811     FORMAT(1x,'ELEBD: ERROR AT POINT:',    i10                   ,/,
     &         1x,'       THE BOUNDARY POINT:',i10                   ,/,
     &         1x,'       DOES NOT BELONG TO ELEMENT : ',i10         ,/,
     &         1x,'       POSSIBLE REASONS:'                         ,/,
     &         1x,'       THE BOUNDARY CONDITION FILE IS NOT      '  ,/,
     &         1x,'       RELEVANT TO THE GEOMETRY FILE           '  ,/,
     &         1x,'       OR THE MESH HAS A WRONG TOPOLOGY')
        CALL plante(1)
        stop
      ENDIF
!
      ENDDO ! K1
!
!     COMPUTES IKLBOR : LIKE IKLE FOR BOUNDARY POINTS, WITH BOUNDARY
!                       POINTS NUMBERING
!
!                       IN 3D WILL BE REDONE IN OTHER ELEB.. ROUTINES
!
!     IF(IELM.NE.41.AND.IELM.NE.51) THEN
      DO k=1,nptfr
        iklbor(k,1) = k
        iklbor(k,2) = kp1bor(k,1)
      ENDDO
!     ENDIF
!
!-----------------------------------------------------------------------
!
!     VERSION 7.0: NOW RENUMBERING IN PARALLEL TO AVOID HOLES IN SEGMENT
!                  NUMBERING. NELEB BECOMES THE REAL NUMBER OF BOUNDARY
!                  ELEMENTS
!
      IF(ncsize.GT.1) THEN
        neleb=0
        DO k=1,nptfr
          IF(kp1bor(k,1).NE.k) THEN
            neleb=neleb+1
            nelbor(neleb)  =nelbor(k)
            nulone(neleb,1)=nulone(k,1)
            nulone(neleb,2)=nulone(k,2)
            iklbor(neleb,1)=iklbor(k,1)
            iklbor(neleb,2)=iklbor(k,2)
          ENDIF
        ENDDO
      ELSE
        neleb=nptfr
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END