flux_ef_vf_3d.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\flux_ef_vf_3d.f
00002 !
00079                      SUBROUTINE FLUX_EF_VF_3D
00080 !                    ************************
00081 !
00082      &(FLOW,W2D,W3D,NSEG2D,NSEG3D,NELEM2,NELEM3,MESH2D,INIFLO,
00083      & IOPT,SENS,IELM3,NPLAN,IKLE,NELMAX,KNOLG)
00084 !
00085 !***********************************************************************
00086 ! BIEF   V6P2                                   21/08/2010
00087 !***********************************************************************
00088 !
00089 !
00090 !
00091 !
00092 !
00093 !
00094 !
00095 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00096 !| FLOW           |<--| FLUX
00097 !| IELM3          |-->| DISCRETISATION IN 3D
00098 !|                |   | (41: PRISMS 51:PRISMS CUT INTO TETRAHEDRONS)
00099 !| IKLE           |-->| CONNECTIVITY TABLE
00100 !| INIFLO         |-->| IF(YES) FLOW WILL BE INITIALISED AT 0.
00101 !| IOPT           |-->| CHOICE OF THE CONSTANT IN FLUX_EF_VF
00102 !| KNOLG          |-->| GIVES THE ORIGINAL GLOBAL NUMBER OF POINTS
00103 !|                |   | IN SCALAR MODE (SIZE NPOIN3 BUT FILLED ONLY
00104 !|                |   | UP TO NPOIN2)
00105 !| MESH2D         |-->| 2D MESH
00106 !| NELEM2         |-->| NUMBER OF ELEMENTS IN 2D
00107 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS (SEE IKLE)
00108 !| NELEM3         |-->| NUMBER OF ELEMENTS IN 3D
00109 !| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH
00110 !| NSEG2D         |-->| NUMBER OF SEGMENTS IN 2D
00111 !| NSEG3D         |-->| NUMBER OF SEGMENTS IN 3D
00112 !| SENS           |-->| IF 1: HORIZONTAL FLUXES FROM BOTTOM TO TOP
00113 !|                |   | IF 2: HORIZONTAL FLUXES FROM TOP TO BOTTOM
00114 !| W2D            |<--| NON ASSEMBLED FLUXES LEAVING POINTS,PER TRIANGLE
00115 !| W3D            |-->| NON ASSEMBLED FLUXES LEAVING POINTS,PER PRISM
00116 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00117 !
00118       USE BIEF, EX_FLUX_EF_VF_3D => FLUX_EF_VF_3D
00119       USE DECLARATIONS_TELEMAC, ONLY : TETRA
00120 !
00121       IMPLICIT NONE
00122       INTEGER LNG,LU
00123       COMMON/INFO/LNG,LU
00124 !
00125 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00126 !
00127       INTEGER, INTENT(IN)             :: NSEG2D,NSEG3D,NELEM2,NELEM3
00128 !                                             *=NSEG2D*NPLAN+NPOIN2*NETAGE
00129       INTEGER, INTENT(IN)             :: IOPT,SENS,IELM3,NPLAN,NELMAX
00130 !                                                    6 IF IELM3=41
00131 !                                                    4 IF IELM3=51
00132       INTEGER, INTENT(IN)             :: IKLE(NELMAX,*),KNOLG(*)
00133       DOUBLE PRECISION, INTENT(INOUT) :: FLOW(NSEG3D)
00134 !                                                   6 IF IELM3=41
00135 !                                                   4 IF IELM3=51
00136       DOUBLE PRECISION, INTENT(IN)    :: W3D(NELEM3,*)
00137       DOUBLE PRECISION, INTENT(INOUT) :: W2D(NELEM2,3)
00138       LOGICAL, INTENT(IN)             :: INIFLO
00139       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH2D
00140 !
00141 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00142 !
00143       INTEGER IPLAN,N1,N2,N3,N4,IELEM,I1,I2,I3,S1,S2,S3,IT,IELEM3D,K,L
00144 !
00145 !     TETRA : WILL GIVE THE LOCAL NUMBERS OF POINTS IN THE PRISM
00146 !             THE 0 CORRESPOND TO SITUATIONS
00147 !             THAT NEVER HAPPEN (TETRA(1,1,1,... OR TETRA(2,2,2,...)
00148 !     INTEGER TETRA(2,2,2,3,4)
00149 !     DATA TETRA / 0,1,1,1,1,1,1,0,0,4,4,4,4,4,4,0,0,6,4,5,5,4,6,0,
00150 !    &             0,2,2,2,2,2,2,0,0,6,6,6,6,6,6,0,0,3,1,2,2,1,3,0,
00151 !    &             0,3,3,3,3,3,3,0,0,5,5,5,5,5,5,0,0,2,3,4,1,6,5,0,
00152 !    &             0,4,5,4,6,6,5,0,0,2,3,3,1,2,1,0,0,4,5,3,6,2,1,0 /
00153 !
00154 !-----------------------------------------------------------------------
00155 !
00156       IF(SENS.NE.1.AND.SENS.NE.2) THEN
00157         IF(LNG.EQ.1) WRITE(LU,*) 'FLUX_EF_VF_3D : SENS INCONNU = ',SENS
00158         IF(LNG.EQ.2) WRITE(LU,*) 'FLUX_EF_VF_3D: UNKNOWN SENS = ',SENS
00159         CALL PLANTE(1)
00160         STOP
00161       ENDIF
00162 !
00163       IF(SENS.EQ.2.AND..NOT.INIFLO) THEN
00164         IF(LNG.EQ.1) THEN
00165           WRITE(LU,*) 'FLUX_EF_VF_3D : SENS = 2 ET INIFLO = .FALSE.'
00166           WRITE(LU,*) '                OPTIONS INCOMPATIBLES'
00167         ENDIF
00168         IF(LNG.EQ.2) THEN
00169           WRITE(LU,*) 'FLUX_EF_VF_3D: SENS = 2 AND INIFLO = .FALSE.'
00170           WRITE(LU,*) '               INCOMPATIBLE OPTIONS'
00171         ENDIF
00172         CALL PLANTE(1)
00173         STOP
00174       ENDIF
00175 !
00176 !-----------------------------------------------------------------------
00177 !
00178       IF(IELM3.EQ.41) THEN
00179 !
00180 !     CASE OF PRISMS
00181 !
00182 !     ADDS FLUXES ON HORIZONTAL SEGMENTS (FOR SEGMENT NUMBERING IN 3D
00183 !                                         SEE STOSEG41 IN BIEF)
00184 !
00185 !     FOR INTERMEDIATE PLANES
00186 !     HERE WE HAVE THE CHOICE OF COMPUTING THE FLUXES BEFORE OR AFTER
00187 !     ASSEMBLING ON VERTICAL. UP TO VERSION 6.1, ASSEMBLY WAS DONE FIRST
00188 !     FROM VERSION 6.2 ON WE COMPUTE THE FLUXES FIRST AND ASSEMBLE AFTER
00189 !     THIS IS MORE COMPATIBLE WITH DISTRIBUTIVE SCHEMES THAT COMPUTE
00190 !     THEIR FLUXES AT ELEMENT LEVEL. THEN THE FLUX LIMITATION DONE
00191 !     IN NA_FLUX3D_LIM WILL BE COMPATIBLE.
00192 !
00193       DO IPLAN=1,NPLAN
00194 !
00195         IF(SENS.EQ.1) THEN
00196           N1=(IPLAN-1    )*NSEG2D+1
00197           N2= IPLAN       *NSEG2D
00198         ELSE
00199           N1=(  NPLAN-IPLAN)*NSEG2D+1
00200           N2=(1+NPLAN-IPLAN)*NSEG2D
00201         ENDIF
00202 !
00203 !       POINTS 1, 2 AND 3 OF UPPER LEVEL
00204         IF(IPLAN.EQ.1) THEN
00205 !         FIRST PLANE: ONLY POINTS 1, 2 AND 3 OF UPPER LEVEL
00206           DO IELEM=1,NELEM2
00207             W2D(IELEM,1)=W3D(IELEM,1)
00208             W2D(IELEM,2)=W3D(IELEM,2)
00209             W2D(IELEM,3)=W3D(IELEM,3)
00210           ENDDO
00211           CALL FLUX_EF_VF(FLOW(N1:N2),W2D,NSEG2D,NELEM2,
00212      &                    MESH2D%ELTSEG%I,MESH2D%ORISEG%I,
00213      &                    MESH2D%IKLE%I,INIFLO,IOPT)
00214         ELSEIF(IPLAN.EQ.NPLAN) THEN
00215 !         LAST PLANE: ONLY POINTS 4, 5 AND 6 OF LOWER LEVEL
00216           N3=NELEM2*(IPLAN-2)
00217           DO IELEM=1,NELEM2
00218             W2D(IELEM,1)=W3D(N3+IELEM,4)
00219             W2D(IELEM,2)=W3D(N3+IELEM,5)
00220             W2D(IELEM,3)=W3D(N3+IELEM,6)
00221           ENDDO
00222           CALL FLUX_EF_VF(FLOW(N1:N2),W2D,NSEG2D,NELEM2,
00223      &                    MESH2D%ELTSEG%I,MESH2D%ORISEG%I,
00224      &                    MESH2D%IKLE%I,INIFLO,IOPT)
00225         ELSE
00226 !       INTERMEDIATE PLANE
00227 !         POINTS 4, 5 AND 6 OF LOWER LEVEL
00228           N3=NELEM2*(IPLAN-2)
00229           DO IELEM=1,NELEM2
00230             W2D(IELEM,1)=W3D(N3+IELEM,4)
00231             W2D(IELEM,2)=W3D(N3+IELEM,5)
00232             W2D(IELEM,3)=W3D(N3+IELEM,6)
00233           ENDDO
00234           CALL FLUX_EF_VF(FLOW(N1:N2),W2D,NSEG2D,NELEM2,
00235      &                    MESH2D%ELTSEG%I,MESH2D%ORISEG%I,
00236      &                    MESH2D%IKLE%I,INIFLO,IOPT)
00237 !         POINTS 1, 2, 3 OF UPPER LEVEL
00238           N4=N3+NELEM2
00239           DO IELEM=1,NELEM2
00240             W2D(IELEM,1)=W3D(N4+IELEM,1)
00241             W2D(IELEM,2)=W3D(N4+IELEM,2)
00242             W2D(IELEM,3)=W3D(N4+IELEM,3)
00243           ENDDO
00244           CALL FLUX_EF_VF(FLOW(N1:N2),W2D,NSEG2D,NELEM2,
00245      &                    MESH2D%ELTSEG%I,MESH2D%ORISEG%I,
00246      &                    MESH2D%IKLE%I,.FALSE.,IOPT)
00247 !                                       !!!!!!!
00248         ENDIF
00249 !
00250       ENDDO
00251 !
00252 !-----------------------------------------------------------------------
00253 !
00254       ELSEIF(IELM3.EQ.51) THEN
00255 !
00256 !     CASE OF PRISMS CUT INTO TETRAHEDRONS
00257 !
00258 !     ADDS FLUXES ON HORIZONTAL SEGMENTS (FOR SEGMENT NUMBERING IN 3D
00259 !                                         SEE STOSEG51 IN BIEF)
00260 !
00261 !
00262 !     INITIALISING
00263 !
00264       DO IELEM=1,NELEM2
00265         W2D(IELEM,1)=0.D0
00266         W2D(IELEM,2)=0.D0
00267         W2D(IELEM,3)=0.D0
00268       ENDDO
00269 !
00270 !     LOOP ON PLANES
00271 !
00272       DO IPLAN=1,NPLAN
00273 !
00274 !       LOOP ON TRIANGLES
00275 !
00276         DO IELEM=1,NELEM2
00277 !
00278 !         HERE LOWER LEVEL OF ELEMENTS ALWAYS TAKEN
00279 !         TO FIND THE WAY THE PRISM HAS BEEN CUT BY LOOKING
00280 !         AT GLOBAL NUMBERS OF POINTS
00281 !         THIS PART IS THUS COMMON TO ALL PLANES
00282 !         IKLE 3D IS TAKEN, COULD BE IKLE 2D AS WELL
00283 !         SO NO DEPENDENCE ON LEVEL IPLAN
00284           I1=IKLE(IELEM,1)
00285           I2=IKLE(IELEM,2)
00286           I3=IKLE(IELEM,3)
00287 !         IN PARALLELISM BACK TO ORIGINAL GLOBAL NUMBERS
00288           IF(NCSIZE.GT.1) THEN
00289             I1=KNOLG(I1)
00290             I2=KNOLG(I2)
00291             I3=KNOLG(I3)
00292           ENDIF
00293 !         THIS IS DONE LIKE IN CPIKLE3 TO USE ARRAY TETRA
00294           IF(I1.GT.I2) THEN
00295             S1=1
00296           ELSE
00297             S1=2
00298           ENDIF
00299           IF(I2.GT.I3) THEN
00300             S2=1
00301           ELSE
00302             S2=2
00303           ENDIF
00304           IF(I3.GT.I1) THEN
00305             S3=1
00306           ELSE
00307             S3=2
00308           ENDIF
00309 !
00310 !         NOW TAKING CONTRIBUTIONS OF TETRAHEDRON K= 1, 2 AND 3
00311 !         OF LOWER AND UPPER LEVEL
00312 !
00313 !         LOWER LEVEL : POINTS 4,5,6 OF ORIGINAL PRISM CONTRIBUTE
00314 !
00315           IF(IPLAN.GT.1) THEN
00316 !           LOOP ON 3 TETRAHEDRONS
00317             DO K=1,3
00318               IELEM3D=3*(IPLAN-2)*NELEM2+IELEM+(K-1)*NELEM2
00319 !             POINTS 1 TO 4
00320               DO L=1,4
00321 !               ORIGINAL NUMBER IN THE PRISM
00322                 IT=TETRA(S1,S2,S3,K,L)
00323                 IF(IT.GE.4) THEN
00324                   W2D(IELEM,IT-3)=W2D(IELEM,IT-3)+W3D(IELEM3D,L)
00325                 ENDIF
00326               ENDDO
00327             ENDDO
00328           ENDIF
00329 !
00330 !         UPPER LEVEL : POINTS 1,2,3 OF ORIGINAL PRISM CONTRIBUTE
00331 !
00332           IF(IPLAN.LT.NPLAN) THEN
00333 !           LOOP ON 3 TETRAHEDRONS
00334             DO K=1,3
00335               IELEM3D=3*(IPLAN-1)*NELEM2+IELEM+(K-1)*NELEM2
00336 !             POINTS 1 TO 4
00337               DO L=1,4
00338 !               ORIGINAL NUMBER IN THE PRISM
00339                 IT=TETRA(S1,S2,S3,K,L)
00340                 IF(IT.LE.3) THEN
00341                   W2D(IELEM,IT)=W2D(IELEM,IT)+W3D(IELEM3D,L)
00342                 ENDIF
00343               ENDDO
00344             ENDDO
00345           ENDIF
00346 !
00347         ENDDO
00348 !
00349         IF(SENS.EQ.1) THEN
00350           N1=(IPLAN-1    )*NSEG2D+1
00351           N2= IPLAN       *NSEG2D
00352         ELSE
00353           N1=(  NPLAN-IPLAN)*NSEG2D+1
00354           N2=(1+NPLAN-IPLAN)*NSEG2D
00355         ENDIF
00356 !
00357       ENDDO
00358 !
00359 !     HERE 2D FLUXES !!!!!!!!!!!!!
00360 !
00361       CALL FLUX_EF_VF(FLOW,W2D,NSEG2D,NELEM2,
00362      &                MESH2D%ELTSEG%I,MESH2D%ORISEG%I,
00363      &                MESH2D%IKLE%I,INIFLO,IOPT)
00364 !
00365 !-----------------------------------------------------------------------
00366 !
00367       ELSE
00368         WRITE(LU,*) 'UNEXPECTED ELEMENT TYPE IN FLUX_EF_VF_3D'
00369         CALL PLANTE(1)
00370         STOP
00371       ENDIF
00372 !
00373 !-----------------------------------------------------------------------
00374 !
00375       RETURN
00376       END

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