coupe.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\coupe.f
00002 !
00059                      SUBROUTINE COUPE
00060 !                    ****************
00061 !
00062      & (F,FINT,SCURV,NPOIN,IKLE3,IFABOR,X,Y,Z,SURFAC,NELEM2,NPOIN3,
00063      &  NETAGE,X1,Y1,Z1,X2,Y2,Z2)
00064 !
00065 !***********************************************************************
00066 ! TELEMAC3D   V6P1                                   21/08/2010
00067 !***********************************************************************
00068 !
00069 !
00070 !
00071 !
00072 !
00073 !
00074 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00075 !| F              |-->| FUNCTION TO INTERPOLATE
00076 !| FINT           |<->| INTERPOLATED FUNCTION ALONG THE PROFILE
00077 !| IFABOR         |-->| CORRESPONDENCE BETWEEN BOUNDARY EDGE AND 2D ELEMENT
00078 !| IKLE3          |-->| GLOBAL 3D CONNECTIVITY
00079 !| NELEM2         |-->| NUMBER OF ELEMENTS IN 2D
00080 !| NETAGE         |-->| NUMBER OF PLANES - 1
00081 !| NPOIN          |-->| NUMBER OF POINTS ALONG THE PROFILE
00082 !| NPOIN3         |-->| NUMBER OF 3D POINTS
00083 !| SCURV          |<->| ABCISSAE OF INTERPOLATING POINTS ON THE PROFILE
00084 !| SURFAC         |-->| AREA OF TRIANGLES IN 2D
00085 !| X              |-->| MESH COORDINATE
00086 !| X1             |-->| COORDINATE OF THE VERTICES OF THE SEGMENT
00087 !| X2             |-->| COORDINATE OF THE VERTICES OF THE SEGMENT
00088 !| Y              |-->| MESH COORDINATE
00089 !| Y1             |-->| COORDINATE OF THE VERTICES OF THE SEGMENT
00090 !| Y2             |-->| COORDINATE OF THE VERTICES OF THE SEGMENT
00091 !| Z              |-->| MESH COORDINATE
00092 !| Z1             |-->| COORDINATE OF THE VERTICES OF THE SEGMENT
00093 !| Z2             |-->| COORDINATE OF THE VERTICES OF THE SEGMENT
00094 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00095 !
00096       IMPLICIT NONE
00097       INTEGER LNG,LU
00098       COMMON/INFO/LNG,LU
00099 !
00100 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00101 !
00102       INTEGER, INTENT(IN) ::  NPOIN, NELEM2, NPOIN3, NETAGE
00103 !
00104       INTEGER, INTENT(IN) ::  IKLE3(NELEM2,NETAGE,6)
00105       INTEGER, INTENT(IN) ::  IFABOR(NELEM2,3)
00106 !
00107       DOUBLE PRECISION, INTENT(IN)    :: F(NPOIN3)
00108       DOUBLE PRECISION, INTENT(INOUT) :: FINT(NPOIN)
00109       DOUBLE PRECISION, INTENT(INOUT) :: SCURV(NPOIN)
00110       DOUBLE PRECISION, INTENT(IN)    :: X(NPOIN3), Y(NPOIN3), Z(NPOIN3)
00111       DOUBLE PRECISION, INTENT(IN)    :: SURFAC(NELEM2)
00112       DOUBLE PRECISION, INTENT(IN)    :: X1,Y1,Z1,X2,Y2,Z2
00113 !
00114 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00115 !
00116 ! LOGICAL VALUE WHICH INDICATES WHETHER THE SUPPORT LINE IS IN THE DOMAIN
00117 !
00118       LOGICAL ENTREE
00119 !
00120       INTEGER IETAGE,IELAR,IPOIN,IFAC,IELARN, NOEUD
00121 !
00122       DOUBLE PRECISION XET,YET,ZET,EPSILO
00123       DOUBLE PRECISION XC,YC,ZC,XS1,YS1,XS2,YS2,A1,A2,DX2
00124       DOUBLE PRECISION XAR1,XAR2,XAR3,YAR1,YAR2,YAR3
00125       DOUBLE PRECISION ZAR1,ZAR2,ZAR3,ZAR4,ZAR5,ZAR6
00126       DOUBLE PRECISION ZS1H,ZS2H,ZS3H,ZS1B,ZS2B,ZS3B,ZHAU,ZBAS
00127       DOUBLE PRECISION ALPHA,BETA,GAMMA
00128       DOUBLE PRECISION SHP1,SHP2,SHP3,SHP4,SHP5,SHP6
00129       DOUBLE PRECISION F1,F2,F3,F4,F5,F6
00130 !
00131       INTRINSIC MIN,ABS,SQRT
00132 !
00133       DATA EPSILO /1.D-6/
00134 !
00135 !***********************************************************************
00136 !
00137 ! COORDINATES OF THE ENTRY POINT OF THE SUPPORT LINE IN THE DOMAIN
00138 !
00139       XET = X1
00140       YET = Y1
00141       ZET = Z1
00142 !
00143       ENTREE = .FALSE.
00144 !
00145       IELAR = 1
00146 !
00147       CALL OV( 'X=C     ' , SCURV, Y , Z , 1.D+60 , NPOIN )
00148       CALL OV( 'X=C     ' , FINT , Y , Z , 1.D+60 , NPOIN )
00149 !
00150       POINTBOUCLE: DO IPOIN = 1,NPOIN
00151 !
00152 !       COORDINATES OF THE CURRENT POINT
00153 !
00154 100     CONTINUE   ! A GOTO TARGET...
00155 !
00156         XC = XET + (X2-X1)*(IPOIN-1)/(NPOIN-1)
00157         YC = YET + (Y2-Y1)*(IPOIN-1)/(NPOIN-1)
00158         ZC = ZET + (Z2-Z1)*(IPOIN-1)/(NPOIN-1)
00159 !
00160 !       MESH ELEMENT THE CURRENT POINT BELONGS TO
00161 !
00162 !       IDENTIFIES THE 2D ELEMENT FIRST, THEN THE LAYER
00163 !
00164 !       2D LOCATION
00165 !
00166 !       CARTESIAN EQUATION OF A FACE: A1*(X-XS1) + A2*(Y-YS1) = 0
00167 !       WHERE XS1,YS1 ARE THE COORDINATES OF A POINT OF THE FACE
00168 !
00169 !       BELONGS TO ELEMENT IELAR: THE ARRIVAL POINT IS ON THE CORRECT
00170 !       SIDE OF THE THREE FACES OF THE ELEMENT
00171 !
00172 30      CONTINUE   ! A GOTO TARGET...
00173 !
00174 !       COORDINATES OF THE FIRST VERTEX
00175         XS1 = X    (IKLE3(IELAR,1     ,1))
00176         YS1 = Y    (IKLE3(IELAR,1     ,1))
00177 !
00178         FACBOUCLE: DO IFAC = 3 , 1 , -1
00179 !
00180 !         ASSUMES VERTICES ARE GIVEN ANTI-CLOCKWISE :
00181 !         USES THE FACT THAT POINT 1 OF FACE N IS
00182 !         POINT 2 OF FACE N-1
00183 !
00184 !         FACE THROUGH VERTICES S1 AND S2
00185           XS2 = XS1
00186           YS2 = YS1
00187           XS1 = X    (IKLE3(IELAR,1     ,IFAC  ))
00188           YS1 = Y    (IKLE3(IELAR,1     ,IFAC  ))
00189 !
00190           A1 = YS1 - YS2
00191           A2 = XS2 - XS1
00192           DX2 = A1*(XC-XS1) + A2*(YC-YS1)
00193 !
00194 !         TEST: IS THE ARRIVAL POINT ON THE CORRECT SIDE OF THE FACE?
00195           IF (DX2.LT.-SURFAC(IELAR)*EPSILO) THEN
00196 !           ON THE WRONG SIDE OF THE FACE
00197 !           GO TO NEXT ELEMENT
00198             IELARN = IFABOR(IELAR,IFAC)
00199             IF (IELARN.LE.0) THEN
00200 !
00201 !             OUTSIDE OF THE DOMAIN
00202               IF (.NOT.ENTREE) THEN
00203 !               IF HAS NOT ALREADY ENTERED, STARTS AGAIN
00204 !               PROGRESSING ALONG THE SUPPORT LINE
00205                 XET = XET + (X2-X1)/(NPOIN-1)
00206                 YET = YET + (Y2-Y1)/(NPOIN-1)
00207                 ZET = ZET + (Z2-Z1)/(NPOIN-1)
00208                 IF ((XET-X2)*(XET-X1).GT.EPSILO*SURFAC(IELAR) .OR.
00209      &             (YET-Y2)*(YET-Y1).GT.EPSILO*SURFAC(IELAR) .OR.
00210      &             (ZET-Z2)*(ZET-Z1).GT.EPSILO*SURFAC(IELAR)) THEN
00211 !                 THE SUPPORT LINE DOES NOT CROSS THE DOMAIN
00212                   IF (LNG.EQ.1) WRITE(LU,21) X1,Y1,Z1,X2,Y2,Z2
00213                   IF (LNG.EQ.2) WRITE(LU,22) X1,Y1,Z1,X2,Y2,Z2
00214                   RETURN
00215                 ENDIF
00216                 GOTO 100
00217               ELSE
00218 !               HAS JUST LEFT THE DOMAIN
00219                 DO NOEUD=IPOIN,NPOIN
00220                   SCURV(NOEUD) = SCURV(IPOIN-1)
00221                   FINT (NOEUD) = FINT (IPOIN-1)
00222                 END DO
00223                 RETURN
00224               ENDIF
00225 !
00226             ELSE
00227               IELAR = IELARN
00228             ENDIF
00229             GOTO 30
00230           ENDIF
00231 !
00232         END DO FACBOUCLE
00233 !
00234 !       IDENTIFIES THE LAYER
00235 !
00236 !       COORDINATES IN THE 2D REFERENCE ELEMENT
00237         XAR1 = X    (IKLE3(IELAR,1     ,1))
00238         XAR2 = X    (IKLE3(IELAR,1     ,2))
00239         XAR3 = X    (IKLE3(IELAR,1     ,3))
00240         YAR1 = Y    (IKLE3(IELAR,1     ,1))
00241         YAR2 = Y    (IKLE3(IELAR,1     ,2))
00242         YAR3 = Y    (IKLE3(IELAR,1     ,3))
00243         ALPHA = ((XC  -XAR1)*(YAR3-YAR1) - (XAR3-XAR1)*(YC  -YAR1))/
00244      &          ((XAR2-XAR1)*(YAR3-YAR1) - (XAR3-XAR1)*(YAR2-YAR1))
00245         BETA  = ((XAR2-XAR1)*(YC  -YAR1) - (XC  -XAR1)*(YAR2-YAR1))/
00246      &          ((XAR2-XAR1)*(YAR3-YAR1) - (XAR3-XAR1)*(YAR2-YAR1))
00247 !
00248 !       COMPUTES THE LAYER
00249         ZS1H = Z(IKLE3(IELAR,NETAGE,4))
00250         ZS2H = Z(IKLE3(IELAR,NETAGE,5))
00251         ZS3H = Z(IKLE3(IELAR,NETAGE,6))
00252         ZS1B = Z(IKLE3(IELAR,1     ,1))
00253         ZS2B = Z(IKLE3(IELAR,1     ,2))
00254         ZS3B = Z(IKLE3(IELAR,1     ,3))
00255 !
00256         ZHAU = (1.D0-ALPHA-BETA)*ZS1H + ALPHA*ZS2H + BETA*ZS3H
00257         ZBAS = (1.D0-ALPHA-BETA)*ZS1B + ALPHA*ZS2B + BETA*ZS3B
00258 !
00259 !       ARE WE IN THE DOMAIN ?
00260 !
00261         IF ( (ZC-ZHAU)/(ZHAU-ZBAS).GT.EPSILO .OR.
00262      &      (ZBAS-ZC)/(ZHAU-ZBAS).GT.EPSILO      ) THEN
00263           IF (.NOT.ENTREE) THEN
00264 !           IF HAS NOT ALREADY ENTERED, STARTS AGAIN
00265 !           PROGRESSING ALONG THE SUPPORT LINE
00266             XET = XET + (X2-X1)/(NPOIN-1)
00267             YET = YET + (Y2-Y1)/(NPOIN-1)
00268             ZET = ZET + (Z2-Z1)/(NPOIN-1)
00269             IF ((XET-X2)*(XET-X1).GT.EPSILO*SURFAC(IELAR) .OR.
00270      &         (YET-Y2)*(YET-Y1).GT.EPSILO*SURFAC(IELAR) .OR.
00271      &         (ZET-Z2)*(ZET-Z1).GT.EPSILO*SURFAC(IELAR)     ) THEN
00272 !             THE SUPPORT LINE DOES NOT CROSS THE DOMAIN
00273               IF (LNG.EQ.1) WRITE(LU,21) X1,Y1,Z1,X2,Y2,Z2
00274               IF (LNG.EQ.2) WRITE(LU,22) X1,Y1,Z1,X2,Y2,Z2
00275               RETURN
00276             ENDIF
00277             GOTO 100
00278           ELSE
00279 !           HAS JUST LEFT THE DOMAIN
00280             DO NOEUD=IPOIN,NPOIN
00281               SCURV(NOEUD) = SCURV(IPOIN-1)
00282               FINT (NOEUD) = FINT (IPOIN-1)
00283             END DO
00284             RETURN
00285           ENDIF
00286         ENDIF
00287 !
00288 !       NUMBER OF THE LAYER
00289         IETAGE = MIN(1+INT(NETAGE*ABS((ZC-ZBAS)/(ZHAU-ZBAS))),NETAGE)
00290 !
00291 !       WE ARE IN THE DOMAIN
00292         IF (.NOT.ENTREE) ENTREE = .TRUE.
00293 !
00294         ZAR1 = Z    (IKLE3(IELAR,IETAGE,1))
00295         ZAR2 = Z    (IKLE3(IELAR,IETAGE,2))
00296         ZAR3 = Z    (IKLE3(IELAR,IETAGE,3))
00297         ZAR4 = Z    (IKLE3(IELAR,IETAGE,4))
00298         ZAR5 = Z    (IKLE3(IELAR,IETAGE,5))
00299         ZAR6 = Z    (IKLE3(IELAR,IETAGE,6))
00300 !
00301 !       COORDINATES IN THE 3D REFERENCE ELEMENT
00302         A1 = ((1.D0-ALPHA-BETA)*ZAR1 + ALPHA*ZAR2 + BETA*ZAR3)/2.D0
00303         A2 = ((1.D0-ALPHA-BETA)*ZAR4 + ALPHA*ZAR5 + BETA*ZAR6)/2.D0
00304         GAMMA = ( ZC-A1-A2)/(A2-A1)
00305 !
00306 !       BASE FUNCTIONS IN THE REFERENCE ELEMENT
00307         SHP1 = (1.D0-ALPHA-BETA)*(1.D0-GAMMA)/2.D0
00308         SHP2 =       ALPHA      *(1.D0-GAMMA)/2.D0
00309         SHP3 =             BETA *(1.D0-GAMMA)/2.D0
00310         SHP4 = (1.D0-ALPHA-BETA)*(1.D0+GAMMA)/2.D0
00311         SHP5 =       ALPHA      *(1.D0+GAMMA)/2.D0
00312         SHP6 =             BETA *(1.D0+GAMMA)/2.D0
00313 !
00314 !       INTERPOLATES
00315         F1 = F    (IKLE3(IELAR,IETAGE,1))
00316         F2 = F    (IKLE3(IELAR,IETAGE,2))
00317         F3 = F    (IKLE3(IELAR,IETAGE,3))
00318         F4 = F    (IKLE3(IELAR,IETAGE,4))
00319         F5 = F    (IKLE3(IELAR,IETAGE,5))
00320         F6 = F    (IKLE3(IELAR,IETAGE,6))
00321 !
00322         FINT(IPOIN) = SHP1*F1 + SHP2*F2 + SHP3*F3 + SHP4*F4 + SHP5*F5
00323      &              + SHP6*F6
00324 !
00325 !       COMPUTES THE CURVILINEAR X-COORDINATE
00326         SCURV(IPOIN) = SQRT( (XC-X1)**2 + (YC-Y1)**2 + (ZC-Z1)**2 )
00327 !
00328       END DO POINTBOUCLE
00329 !
00330 !-----------------------------------------------------------------------
00331 !
00332 21    FORMAT('COUPE : LA COUPE NE TRAVERSE PAS LE DOMAINE',/,
00333      &       'X1,Y1,Z1 :',3G16.7,/,'X2,Y2,Z2 :',3G16.7)
00334 22    FORMAT('COUPE : THE CROSSING LINE DO NOT CROSS THE DOMAIN',/,
00335      &       'X1,Y1,Z1 :',3G16.7,/,'X2,Y2,Z2 :',3G16.7)
00336 !
00337 !-----------------------------------------------------------------------
00338 !
00339       RETURN
00340       END SUBROUTINE COUPE

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