infcel.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\infcel.f
00002 !
00061                      SUBROUTINE INFCEL
00062 !                    *****************
00063 !
00064      &(XX,YY,NUBO,VNOIN,NPOIN,NELEM,NSEG,CMI,AIRST,
00065      & GLOSEG,COORD_G,ELTSEG,ORISEG,IFABOR)
00066 !
00067 !***********************************************************************
00068 ! BIEF   V6P3                                   21/12/20112
00069 !***********************************************************************
00070 !
00071 !          COMPUTES  CMI, AIRST AND NORMALS (VNOIN)
00072 !
00073 !
00074 !
00075 !
00076 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00077 !| IFAPAR         |-->| IFAPAR(1:3,IELEM)=PROCESSOR NUMBERS BEHIND THE
00078 !|                |   | 3 ELEMENT EDGES  (NUMBERS FROM 0 TO NCSIZE-1)
00079 !|                |   | IFAPAR(4:6,IELEM): -LOCAL- ELEMENT NUMBERS
00080 !|                |   | BEHIND THE 3 EDGES
00081 !| AIRST          |<--| AREAS OF SUB-TRIANGLES IG1G2
00082 !| CMI            |<--| COORDINATES OF MID-INTERFACE POINTS
00083 !| INDPU          |-->| INDEX TABLE : IF 0: NOT AN INTERFACE POINT
00084 !|                |   |               IF NOT 0: ADDRESS IN THE LIST
00085 !|                |   |               OF BOUNDARY POINTS.
00086 !| IFABOR         |-->| IFABOR(IEL,I) IS THE ELEMENT BEHIND THE EDGE I OF
00087 !|                |---|  ELEMENT IEL, OTHERWISE
00088 !|                |---|     IFABOR(IEL,I) = -2 : THIS IS INTERFACE EDGE
00089 !|                |---|     IFABOR(IEL,I) = 0  : THIS IS BOUNDARY EDGE
00090 !|                |---|     IFABOR(IEL,I) = -1 : THIS IS LIQUID BOUNDARY EDGE
00091 !| KNOLG          |-->| CONVERSION FROM LOCAL TO GLOBAL NUMEBEING
00092 !| KP1BOR         |<--| NUMBER OF FOLLOWING AND PRECEDING NODE OF BOUNDARY POINT K.
00093 !| MESH           |-- | MESH STRUCTURE
00094 !| NB_NEIGHB_SEG  |-->| NUMBER OF NEIGHBOURING SUB-DOMAINS (FOR EDGES)
00095 !| NELEM          |-->| NUMBER OF ELEMENTS
00096 !| NPOIN          |-->| NUMBER OF POINTS
00097 !| NSEG           |-->| NUMBER OF SEGMENTS
00098 !| NUBO           |<--| FIRST AND SECOND POINT OF SEGMENTS (GLOSEG ?)
00099 !| VNOIN          |<--| NORMAL TO THE INTERFACE
00100 !|                |   | (2 FIRST COMPONENTS) AND
00101 !|                |   | SEGMENT LENGTH (3RD COMPONENT)
00102 !| XX             |-->| ABSCISSAE OF POINTS IN THE MESH
00103 !| YY             |-->| ORDINATES OF POINTS IN THE MESH
00104 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00105 !
00106       USE BIEF_DEF
00107       IMPLICIT  NONE
00108 !
00109       INTEGER LNG,LU
00110       COMMON/INFO/LNG,LU
00111 !
00112 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00113 !
00114       INTEGER, INTENT(IN)             :: NSEG,NPOIN,NELEM
00115 !     TO REPLACE IN NEXT RELEASE
00116       INTEGER, INTENT(INOUT)          :: NUBO(2,NSEG)
00117       INTEGER, INTENT(IN)             :: GLOSEG(NSEG,2)
00118       INTEGER, INTENT(IN)             :: ELTSEG(NELEM,3)
00119       INTEGER, INTENT(IN)             :: ORISEG(NELEM,3)
00120       DOUBLE PRECISION, INTENT(IN)    :: XX(NPOIN),YY(NPOIN),CMI(2,NSEG)
00121       DOUBLE PRECISION, INTENT(INOUT) :: VNOIN(3,NSEG)
00122       DOUBLE PRECISION, INTENT(INOUT) :: AIRST(2,*)
00123       DOUBLE PRECISION, INTENT(IN)    :: COORD_G(NSEG,4)
00124       INTEGER, INTENT(IN)             :: IFABOR(NELEM,3)
00125 !
00126 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00127 !
00128       INTEGER NB1,NB2,ISEG,IEL,FACT,I,IER
00129       DOUBLE PRECISION,PARAMETER :: EPS= 1.D-14
00130       DOUBLE PRECISION X1,Y1,X2,Y2,RNORM,XGG,YGG,XG1,YG1,XG2,YG2
00131 !
00132       INTRINSIC ABS,SQRT
00133 !
00134       LOGICAL, ALLOCATABLE :: YESNO(:)
00135       ALLOCATE(YESNO(NSEG),STAT=IER)
00136       IF(IER.NE.0) THEN
00137         IF(LNG.EQ.1) WRITE(LU,*) 'INFCEL: ERREUR D''ALLOCATION DE YESNO'
00138         IF(LNG.EQ.2) WRITE(LU,*) 'INFCEL: ALLOCATION ERROR OF YESNO'
00139         CALL PLANTE(1)
00140         STOP
00141       ENDIF
00142 !
00143 !-----------------------------------------------------------------------
00144 !
00145 !-----------------
00146 ! 1. INITIALISES
00147 !-----------------
00148 !
00149       DO ISEG = 1 , NSEG
00150         VNOIN(1,ISEG) = 0.D0
00151         VNOIN(2,ISEG) = 0.D0
00152         VNOIN(3,ISEG) = 0.D0
00153         YESNO(ISEG)   =.FALSE.
00154       ENDDO
00155 !
00156 !-----------------
00157 ! 2. MAIN LOOP
00158 !-----------------
00159 !
00160       DO IEL=1, NELEM
00161         DO I = 1,3
00162           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00163             FACT=1
00164             ISEG = ELTSEG(IEL,I)
00165 !           LET S RECUPERATE NODE NUMBERS
00166             NB1 = GLOSEG(ISEG,1)
00167             NB2 = GLOSEG(ISEG,2)
00168 !           TO REMOVE REDUNDANT
00169             NUBO(1,ISEG) = NB1
00170             NUBO(2,ISEG) = NB2
00171 !           THEIR COORDINATES
00172             X1 = XX(NB1)
00173             Y1 = YY(NB1)
00174             X2 = XX(NB2)
00175             Y2 = YY(NB2)
00176 !           CENTER OF GRAVITY OF NEIGHBORING ELEMENTS
00177             XG1 = COORD_G(ISEG,1)
00178             YG1 = COORD_G(ISEG,2)
00179             XG2 = COORD_G(ISEG,3)
00180             YG2 = COORD_G(ISEG,4)
00181             IF(IFABOR(IEL,I).EQ.-1.OR.IFABOR(IEL,I).EQ.0)THEN  ! BOUNDARY SEGMENT
00182               IF(ABS(XG1).LT.EPS.AND.ABS(YG1).LT.EPS)THEN   !BOUNDARY SEGMENT (TO IMPROVE THE TEST!!!)
00183                 XG1=CMI(1,ISEG)                   ! THIS CASE COULD REALLY HAPPEN ?!
00184                 YG1=CMI(2,ISEG)
00185               ELSEIF(ABS(XG2).LT.EPS.AND.ABS(YG2).LT.EPS)THEN
00186                 XG2=CMI(1,ISEG)
00187                 YG2=CMI(2,ISEG)
00188               ENDIF
00189             ENDIF
00190 !           SURFACE OF TRIANGLES (NB1,G1,G2) AND (NB2,G1,G2)
00191             AIRST(1,ISEG)=0.5D0*ABS((X1-XG1)*(Y1-YG2)-
00192      &                              (Y1-YG1)*(X1-XG2))
00193             AIRST(2,ISEG)=0.5D0*ABS((X2-XG1)*(Y2-YG2)-
00194      &                              (Y2-YG1)*(X2-XG2))
00195 !           NORMAL VECTORS TO INTERFACES AND THEIR LENGHT
00196             XGG = XG1-XG2
00197             YGG = YG1-YG2
00198             RNORM=SQRT(XGG**2+YGG**2)
00199             IF(RNORM.GT.EPS) THEN
00200               IF(ORISEG(IEL,I).EQ.2) FACT=-FACT
00201               VNOIN(1,ISEG) =  FACT*YGG/RNORM
00202               VNOIN(2,ISEG) = -FACT*XGG/RNORM
00203               VNOIN(3,ISEG) = RNORM
00204             ELSE
00205               WRITE(LU,*)'**************************************'
00206               IF(LNG.EQ.1) THEN
00207                 WRITE(LU,*)'INFCEL: LONGUEUR D INTERFACE NULLE',RNORM
00208                 WRITE(LU,*)'POUR LE SEGMENT',ISEG
00209                 WRITE(LU,*)'AVEC LES NOEUDS GLOBAUX',NB1,NB2
00210               ELSE
00211                 WRITE(LU,*)'INFCEL: INTERFACE LENGTH NIL',RNORM
00212                 WRITE(LU,*)'FOR SEGMENT',ISEG
00213                 WRITE(LU,*)'WITH GLOBAL NODES',NB1,NB2
00214               ENDIF
00215               WRITE(LU,*)'**************************************'
00216               CALL PLANTE(1)
00217               STOP
00218             ENDIF
00219             YESNO(ISEG)=.TRUE.
00220           ENDIF
00221         ENDDO
00222       ENDDO
00223 !
00224 !---------------------------------------------------------------------
00225 !
00226       RETURN
00227       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0