infcel.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE infcel
!                   *****************
!
     &(xx,yy,nubo,vnoin,npoin,nelem,nseg,cmi,airst,
     & gloseg,coord_g,eltseg,oriseg,ifabor)
!
!***********************************************************************
! BIEF   V6P3                                   21/12/20112
!***********************************************************************
!
!brief     REPLACE  OLD INFCEL: NOW COMMON NUBO (GLOSEG) WITH FE
!          COMPUTES  CMI, AIRST AND NORMALS (VNOIN)
!
!history
!+        18/06/03
!+        V5P4
!+
!
!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 R. ATA (EDF R&D - LNHE)
!+        21/01/2013
!+        V6P3
!+   rewritten for new data structure of finite volumes
!+
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IFAPAR         |-->| IFAPAR(1:3,IELEM)=PROCESSOR NUMBERS BEHIND THE
!|                |   | 3 ELEMENT EDGES  (NUMBERS FROM 0 TO NCSIZE-1)
!|                |   | IFAPAR(4:6,IELEM): -LOCAL- ELEMENT NUMBERS
!|                |   | BEHIND THE 3 EDGES
!| AIRST          |<--| AREAS OF SUB-TRIANGLES IG1G2
!| CMI            |<--| COORDINATES OF MID-INTERFACE POINTS
!| INDPU          |-->| INDEX TABLE : IF 0: NOT AN INTERFACE POINT
!|                |   |               IF NOT 0: ADDRESS IN THE LIST
!|                |   |               OF BOUNDARY POINTS.
!| IFABOR         |-->| IFABOR(IEL,I) IS THE ELEMENT BEHIND THE EDGE I OF
!|                |---|  ELEMENT IEL, OTHERWISE
!|                |---|     IFABOR(IEL,I) = -2 : THIS IS INTERFACE EDGE
!|                |---|     IFABOR(IEL,I) = 0  : THIS IS BOUNDARY EDGE
!|                |---|     IFABOR(IEL,I) = -1 : THIS IS LIQUID BOUNDARY EDGE
!| KNOLG          |-->| CONVERSION FROM LOCAL TO GLOBAL NUMEBEING
!| KP1BOR         |<--| NUMBER OF FOLLOWING AND PRECEDING NODE OF BOUNDARY POINT K.
!| MESH           |-- | MESH STRUCTURE
!| NB_NEIGHB_SEG  |-->| NUMBER OF NEIGHBOURING SUB-DOMAINS (FOR EDGES)
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| NSEG           |-->| NUMBER OF SEGMENTS
!| NUBO           |<--| FIRST AND SECOND POINT OF SEGMENTS (GLOSEG ?)
!| VNOIN          |<--| NORMAL TO THE INTERFACE
!|                |   | (2 FIRST COMPONENTS) AND
!|                |   | SEGMENT LENGTH (3RD COMPONENT)
!| XX             |-->| ABSCISSAE OF POINTS IN THE MESH
!| YY             |-->| ORDINATES OF POINTS IN THE MESH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief_def
      IMPLICIT  NONE
!
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NSEG,NPOIN,NELEM
!     TO REPLACE IN NEXT RELEASE
      INTEGER, INTENT(INOUT)          :: NUBO(2,nseg)
      INTEGER, INTENT(IN)             :: GLOSEG(nseg,2)
      INTEGER, INTENT(IN)             :: ELTSEG(nelem,3)
      INTEGER, INTENT(IN)             :: ORISEG(nelem,3)
      DOUBLE PRECISION, INTENT(IN)    :: XX(npoin),YY(npoin),CMI(2,nseg)
      DOUBLE PRECISION, INTENT(INOUT) :: VNOIN(3,nseg)
      DOUBLE PRECISION, INTENT(INOUT) :: AIRST(2,*)
      DOUBLE PRECISION, INTENT(IN)    :: COORD_G(nseg,4)
      INTEGER, INTENT(IN)             :: IFABOR(nelem,3)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER NB1,NB2,ISEG,IEL,FACT,I,IER
      DOUBLE PRECISION,PARAMETER :: EPS= 1.d-14
      DOUBLE PRECISION X1,Y1,X2,Y2,RNORM,XGG,YGG,XG1,YG1,XG2,YG2
!
      INTRINSIC abs,sqrt
!
      LOGICAL, ALLOCATABLE :: YESNO(:)
      ALLOCATE(yesno(nseg),stat=ier)
      CALL check_allocate(ier, 'INFCEL')
!
!-----------------------------------------------------------------------
!
!-----------------
! 1. INITIALISES
!-----------------
!
      DO iseg = 1 , nseg
        vnoin(1,iseg) = 0.d0
        vnoin(2,iseg) = 0.d0
        vnoin(3,iseg) = 0.d0
        yesno(iseg)   =.false.
      ENDDO
!
!-----------------
! 2. MAIN LOOP
!-----------------
!
      DO iel=1, nelem
        DO i = 1,3
          IF(.NOT.yesno(eltseg(iel,i)))THEN
            fact=1
            iseg = eltseg(iel,i)
!           LET S RECUPERATE NODE NUMBERS
            nb1 = gloseg(iseg,1)
            nb2 = gloseg(iseg,2)
!           TO REMOVE REDUNDANT
            nubo(1,iseg) = nb1
            nubo(2,iseg) = nb2
!           THEIR COORDINATES
            x1 = xx(nb1)
            y1 = yy(nb1)
            x2 = xx(nb2)
            y2 = yy(nb2)
!           CENTER OF GRAVITY OF NEIGHBORING ELEMENTS
            xg1 = coord_g(iseg,1)
            yg1 = coord_g(iseg,2)
            xg2 = coord_g(iseg,3)
            yg2 = coord_g(iseg,4)
            IF(ifabor(iel,i).EQ.-1.OR.ifabor(iel,i).EQ.0)THEN  ! BOUNDARY SEGMENT
              IF(abs(xg1).LT.eps.AND.abs(yg1).LT.eps)THEN   !BOUNDARY SEGMENT (TO IMPROVE THE TEST!!!)
                xg1=cmi(1,iseg)                   ! THIS CASE COULD REALLY HAPPEN ?!
                yg1=cmi(2,iseg)
              ELSEIF(abs(xg2).LT.eps.AND.abs(yg2).LT.eps)THEN
                xg2=cmi(1,iseg)
                yg2=cmi(2,iseg)
              ENDIF
            ENDIF
!           SURFACE OF TRIANGLES (NB1,G1,G2) AND (NB2,G1,G2)
            airst(1,iseg)=0.5d0*abs((x1-xg1)*(y1-yg2)-
     &                              (y1-yg1)*(x1-xg2))
            airst(2,iseg)=0.5d0*abs((x2-xg1)*(y2-yg2)-
     &                              (y2-yg1)*(x2-xg2))
!           NORMAL VECTORS TO INTERFACES AND THEIR LENGHT
            xgg = xg1-xg2
            ygg = yg1-yg2
            rnorm=sqrt(xgg**2+ygg**2)
            IF(rnorm.GT.eps) THEN
              IF(oriseg(iel,i).EQ.2) fact=-fact
              vnoin(1,iseg) =  fact*ygg/rnorm
              vnoin(2,iseg) = -fact*xgg/rnorm
              vnoin(3,iseg) = rnorm
            ELSE
              WRITE(lu,*)'**************************************'
              WRITE(lu,*)'INFCEL: INTERFACE LENGTH NIL',rnorm
              WRITE(lu,*)'FOR SEGMENT',iseg
              WRITE(lu,*)'WITH GLOBAL NODES',nb1,nb2
              WRITE(lu,*)'**************************************'
              CALL plante(1)
              stop
            ENDIF
            yesno(iseg)=.true.
          ENDIF
        ENDDO
      ENDDO
!
!---------------------------------------------------------------------
!
      RETURN
      END