stoseg.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE stoseg
!                   *****************
!
     &(ifabor,nelem,nelmax,nelmax2,ielm,ikle,nbor,nptfr,
     & gloseg,maxseg,eltseg,oriseg,nseg,nelbor,nulone,knolg,
     & iklbor,nelebx,neleb)
!
!***********************************************************************
! BIEF   V7P0                                   21/08/2010
!***********************************************************************
!
!brief    BUILDS THE DATA STRUCTURE FOR EDGE-BASED STORAGE.
!
!history  J-M HERVOUET (LNH)
!+        02/10/08
!+        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 (EDF LAB, LNHE)
!+        13/03/2014
!+        V7P0
!+   Now written to enable different numbering of boundary points and
!+   boundary segments.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ELTSEG         |<--| SEGMENTS OF EVERY TRIANGLE.
!| GLOSEG         |<--| GLOBAL NUMBERS OF POINTS OF SEGMENTS.
!| 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.
!| KNOLG          |-->| GLOBAL NUMBER OF A LOCAL POINT IN PARALLEL
!| MAXSEG         |<--| MAXIMUM NUMBER OF SEGMENTS
!| NBOR           |-->| GLOBAL NUMBERS OF BOUNDARY POINTS.
!| NELBOR         |-->| NUMBER OF ELEMENT CONTAINING SEGMENT K OF
!|                |   | THE BOUNDARY.
!| NELEM          |-->| NUMBER OF ELEMENTS IN THE MESH
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS IN 3D
!| NELMAX2        |-->| MAXIMUM NUMBER OF ELEMENTS IN 2D
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS.
!| NSEG           |<--| NUMBER OF SEGMENTS OF THE MESH.
!| NULONE         |-->| LOCAL NUMBER OF BOUNDARY POINTS IN A BOUNDARY
!|                |   | ELEMENT.
!| ORISEG         |<--| ORIENTATION OF SEGMENTS OF EVERY TRIANGLE.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_stoseg => stoseg
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)    :: NELMAX,NELMAX2,NPTFR,NSEG,MAXSEG,IELM
      INTEGER, INTENT(IN)    :: NELEM,NELEBX,NELEB
      INTEGER, INTENT(IN)    :: NBOR(nptfr)
      INTEGER, INTENT(IN)    :: IFABOR(nelmax2,*),IKLE(nelmax,*)
      INTEGER, INTENT(IN)    :: NELBOR(*),NULONE(*)
      INTEGER, INTENT(IN)    :: IKLBOR(nelebx,2)
      INTEGER, INTENT(INOUT) :: GLOSEG(maxseg,2)
      INTEGER, INTENT(INOUT) :: ELTSEG(nelmax,*),ORISEG(nelmax,3)
      INTEGER, INTENT(IN)    :: KNOLG(*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER NSE,IELEB
!
      INTEGER NEL,IFA,I1,I2,J1,J2,IFACE,JFACE,IG1,IG2
      INTEGER IELEM,IELEM1,IELEM2
!
      INTEGER :: NEXT(3)
      parameter( next = (/ 2,3,1 /) )
!
!-----------------------------------------------------------------------
!
      IF(ielm.NE.11.AND.ielm.NE.12.AND.ielm.NE.13.AND.ielm.NE.14) THEN
        WRITE(lu,501) ielm
501     FORMAT(1x,'STOSEG (BIEF) : UNEXPECTED ELEMENT: ',1i6)
        CALL plante(1)
        stop
      ENDIF
!
!     INITIALISES ELTSEG
!
      DO ielem = 1 , nelem
        eltseg(ielem,1) = 0
        eltseg(ielem,2) = 0
        eltseg(ielem,3) = 0
      ENDDO
!
!-----------------------------------------------------------------------
!
!     LOOP ON BOUNDARY POINTS :
!
      DO ieleb = 1 , neleb
!
!         NOTE: ON BOUNDARIES, SEGMENTS ARE NOT ORIENTED LOWER RANK
!               TO HIGHER RANK, AS IS DONE FOR INTERNAL SEGMENTS
          gloseg(ieleb,1) = nbor(iklbor(ieleb,1))
          gloseg(ieleb,2) = nbor(iklbor(ieleb,2))
          nel = nelbor(ieleb)
          ifa = nulone(ieleb)
          eltseg(nel,ifa) = ieleb
          oriseg(nel,ifa) = 1
!
      ENDDO
!
      nse=neleb
!
!-----------------------------------------------------------------------
!
!     LOOP ON ELEMENTS FOR NUMBERING INTERNAL SEGMENTS AND FILLING:
!     GLOSEG, ELTSEG, ORISEG
!
      DO ielem1 = 1 , nelem
        DO iface = 1 , 3
          IF(eltseg(ielem1,iface).EQ.0) THEN
!           NEW SEGMENT (HENCE INTERNAL SO IFABOR<>0)
            nse = nse + 1
!           BOTH NEIGHBOURING ELEMENTS ARE TREATED FOR THIS SEGMENT
            i1 = ikle(ielem1,     iface)
            i2 = ikle(ielem1,next(iface))
            IF(i1.EQ.i2) THEN
              WRITE(lu,*) 'STOSEG: EDGE MADE OF ONLY ONE POINT'
              WRITE(lu,*) '        ELEMENT ',ielem1,' FACE ',iface
              CALL plante(1)
              stop
            ENDIF
            eltseg(ielem1,iface) = nse
            IF(ncsize.GT.1) THEN
              ig1=knolg(i1)
              ig2=knolg(i2)
            ELSE
              ig1=i1
              ig2=i2
            ENDIF
!           SEGMENT ORIENTED LOWER RANK TO HIGHER RANK
            IF(ig1.LT.ig2) THEN
              gloseg(nse,1) = i1
              gloseg(nse,2) = i2
              oriseg(ielem1,iface) = 1
            ELSE
              gloseg(nse,1) = i2
              gloseg(nse,2) = i1
              oriseg(ielem1,iface) = 2
            ENDIF
!           OTHER ELEMENT NEIGHBOURING THIS SEGMENT
            ielem2 = ifabor(ielem1,iface)
!           IELEM2 = 0 OR -1 MAY OCCUR IN PARALLEL MODE
            IF(ielem2.GT.0) THEN
!             LOOKS FOR THE RIGHT SIDE OF ELEMENT IELEM2
              DO jface = 1,3
                j1 = ikle(ielem2,     jface)
                j2 = ikle(ielem2,next(jface))
!               ALL ELEMENTS HAVE A COUNTER-CLOCKWISE NUMBERING
                IF(i1.EQ.j2.AND.i2.EQ.j1) THEN
                  eltseg(ielem2,jface) = nse
                  oriseg(ielem2,jface) = 3-oriseg(ielem1,iface)
!                 SIDE FOUND, NO NEED TO GO ON
                  GO TO 1000
                ELSEIF(i1.EQ.j1.AND.i2.EQ.j2) THEN
!                 SIDE BADLY ORIENTED
                  WRITE(lu,*) 'STOSEG: WRONG MESH'
                  WRITE(lu,*) '        FACE ',jface
                  WRITE(lu,*) '        OF ELEMENT ',ielem2
                  WRITE(lu,*) '        IS NOT WELL ORIENTED'
                  WRITE(lu,*) '         (POINTS ',i1,' AND ',i2,')'
                  CALL plante(1)
                  stop
                ENDIF
              ENDDO
!             SIDE NOT FOUND, THIS IS AN ERROR
              WRITE(lu,*) 'STOSEG: WRONG MESH'
              WRITE(lu,*) '        ELEMENTS ',ielem1,' AND ',ielem2
              WRITE(lu,*) '        LINKED BY POINTS ',i1,' AND ',i2
              WRITE(lu,*) '        BUT THESE POINTS ARE NOT AN EDGE'
              WRITE(lu,*) '        OF ELEMENT ',ielem2
              CALL plante(1)
              stop
            ENDIF
1000        CONTINUE
          ENDIF
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
!     CHECKS
!
      IF(nseg.NE.nse) THEN
        WRITE(lu,503) nse,nseg
503     FORMAT(1x,'STOSEG (BIEF): WRONG NUMBER OF SEGMENTS : ',1i6,
     &            '               INSTEAD OF ',1i6,' EXPECTED')
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END