mt14tt.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\mt14tt.f
00002 !
00053                      SUBROUTINE MT14TT
00054 !                    *****************
00055 !
00056      &(T,XM,PPQ,LEGO,XMUL,SU,SV,SW,U,V,W,SF,SG,SH,F,G,H,
00057      & X,Y,Z,SURFAC,IKLE,NELEM,NELMAX,SIGMAG,SPECAD,NPLAN,NPOIN2)
00058 !
00059 !***********************************************************************
00060 ! BIEF   V6P2                                   21/08/2010
00061 !***********************************************************************
00062 !
00063 !
00064 !
00065 !reference  J-M HERVOUET THESIS OR BOOK: MURD SCHEME IN 3 DIMENSIONS
00066 !+             CHAPTER 6.
00067 !+
00068 !+          Computational Fluid dynamics '92 volume 1. Multidimensional
00069 !+          upwind schemes for scalar advection on tetrahedral meshes.
00070 !+          G. Bourgois, H. Deconinck, P.L. Roe and R. Struijs
00071 !
00072 !
00073 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00074 !| F              |-->| FUNCTION USED IN THE FORMULA
00075 !| G              |-->| FUNCTION USED IN THE FORMULA
00076 !| H              |-->| FUNCTION USED IN THE FORMULA
00077 !| IKLE           |-->| CONNECTIVITY TABLE
00078 !| LEGO           |-->| LOGICAL. IF YES RESULT ASSEMBLED.
00079 !| NELEM          |-->| NUMBER OF ELEMENTS
00080 !| NELEM2         |-->| NUMBER OF TRIANGLES IN THE UNDERLYING 2D MESH
00081 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00082 !| NPLAN          |-->| NUMBER OF PLANES IN THE MESH
00083 !| NPOIN2         |-->| NUMBER OF POINTS IN THE UNDERLYING 2D MESH
00084 !| PPQ            |-->| STORAGE IN XM OF OFF-DIAGONAL ELEMENTS
00085 !| SF             |-->| BIEF_OBJ STRUCTURE OF F
00086 !| SG             |-->| BIEF_OBJ STRUCTURE OF G
00087 !| SH             |-->| BIEF_OBJ STRUCTURE OF H
00088 !| SIGMAG         |-->| IF YES, GENERALISED SIGMA TRANSFORMATION
00089 !| SPECAD         |-->| IF YES, SPECIAL FORM OF ADVECTION FIELD
00090 !| SU             |-->| BIEF_OBJ STRUCTURE OF U
00091 !| SV             |-->| BIEF_OBJ STRUCTURE OF V
00092 !| SW             |-->| BIEF_OBJ STRUCTURE OF W
00093 !| SURFAC         |-->| AREA OF 2D ELEMENTS
00094 !| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
00095 !| U              |-->| FUNCTION USED IN THE FORMULA
00096 !| V              |-->| FUNCTION USED IN THE FORMULA
00097 !| W              |-->| FUNCTION USED IN THE FORMULA
00098 !| X              |-->| ABSCISSAE OF POINTS IN THE MESH
00099 !| Y              |-->| ORDINATES OF POINTS IN THE MESH
00100 !| Z              |-->| ELEVATIONS OF POINTS IN THE MESH
00101 !| XM             |<->| OFF-DIAGONAL TERMS
00102 !| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
00103 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00104 !
00105       USE BIEF           !, EX_MT14TT => MT14TT
00106       USE DECLARATIONS_TELEMAC, ONLY : ISEGT
00107 !
00108       IMPLICIT NONE
00109       INTEGER LNG,LU
00110       COMMON/INFO/LNG,LU
00111 !
00112 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00113 !
00114       INTEGER, INTENT(IN) :: NELEM,NELMAX,NPLAN,NPOIN2
00115       INTEGER, INTENT(IN) :: IKLE(NELMAX,4),PPQ(6,6)
00116 !
00117       DOUBLE PRECISION, INTENT(INOUT) :: T(NELMAX,4),XM(12,NELMAX)
00118 !
00119       DOUBLE PRECISION, INTENT(IN) :: XMUL
00120       DOUBLE PRECISION, INTENT(IN) :: U(*),V(*),W(*),F(*),G(*)
00121       DOUBLE PRECISION, INTENT(IN) :: H(NELMAX,4)
00122 !
00123       LOGICAL, INTENT(IN) :: LEGO,SIGMAG,SPECAD
00124 !
00125 !     STRUCTURES DE U,V,W
00126 !
00127       TYPE(BIEF_OBJ), INTENT(IN) :: SU,SV,SW,SF,SG,SH
00128 !
00129       DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),Z(*)
00130       DOUBLE PRECISION, INTENT(IN) :: SURFAC(NELMAX)
00131 !
00132 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00133 !
00134 !     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
00135 !
00136       INTEGER IELEM,IPLAN,J1,J2,ISEG,II1,II2,II3,IELEM3D,NELEM2
00137 !
00138       DOUBLE PRECISION X2,X3,Y2,Y3
00139       DOUBLE PRECISION L12,L13,L14,L23,L24,L34,L21,L31,L41,L32,L42,L43
00140       DOUBLE PRECISION ALFA(4),SURFSUR3,SUMMAXK
00141 !
00142       INTEGER IL1,IL2,I1,I2,I3,I4,I5,I6,I
00143 !
00144 !-----------------------------------------------------------------------
00145 !
00146 !     THE SIX SEGMENTS IN A TETRAHEDRON
00147 !     ISEGT(ISEG,1 OR 2) : FIRST OR SECOND POINT OF SEGMENT ISEG
00148 !     INTEGER ISEGT(6,2)
00149 !     DATA ISEGT/1,2,3,1,2,3,2,3,1,4,4,4/
00150 !
00151 !=======================================================================
00152 !
00153 !     RETRIEVING THE NUMBER OF ELEMENTS THE 2D MESH
00154 !
00155       NELEM2=NELMAX/3/(NPLAN-1)
00156 !                                                         *   *
00157 !     COMPUTES THE COEFFICIENTS A(I) (INTEGRAL OF H U.GRAD(PSI), ONLY
00158 !     CONSIDERING THE HORIZONTAL PART OF U). SEE JMH THESIS
00159 !
00160 !     NOTE: THIS VC04TT CALL IS ALREADY DONE BY FLUX3D IN TELEMAC-3D
00161 !           THROUGH A CALL TO VECTOR. THE EQUIVALENT OF T HAS BEEN
00162 !           SAVED IN WHAT IS HERE H
00163 !
00164 !     FORMUL='             HOR'
00165 !     CALL VC04TT(XMUL,SU,SV,SW,U,V,W,SF,SG,SH,F,G,H,X,Y,Z,
00166 !    * IKLE(1,1),IKLE(1,2),IKLE(1,3),IKLE(1,4),
00167 !    * NELEM,NELMAX,T(1,1),T(1,2),T(1,3),T(1,4),
00168 !    * SPECAD,FORMUL,NPLAN)
00169 !
00170 !     NOW T IS RETRIEVED BY STORAGE DONE INTO FLUX3D
00171 !     FUNCTION H MUST CONTAIN THE EQUIVALENT OF T ABOVE, WHICH IS
00172 !     THE NON ASSEMBLED VALUE OF FLUINT
00173 !
00174 !     DO IELEM=1,NELEM
00175 !       T(IELEM,1)=H(         IELEM)
00176 !       T(IELEM,2)=H(  NELMAX+IELEM)
00177 !       T(IELEM,3)=H(2*NELMAX+IELEM)
00178 !       T(IELEM,4)=H(3*NELMAX+IELEM)
00179 !     ENDDO
00180 !
00181 !     FORMUL NORMALLY STARTS WITH VGRADP OR VGRADP2, OR VGRADP22 TO KNOW
00182 !     SIGMAG AND SPECAD, USELESS HERE.
00183 !
00184 !-----------------------------------------------------------------------
00185 !
00186       IF(SW%ELM.NE.51) THEN
00187         WRITE(LU,*) 'MT14TT DISCRETISATION NON PREVUE'
00188         CALL PLANTE(1)
00189         STOP
00190       ENDIF
00191 !
00192 !-----------------------------------------------------------------------
00193 !
00194       DO IPLAN=1,NPLAN-1
00195         DO IELEM=1,NELEM2
00196 !
00197 !         THE SIX GLOBAL NUMBERS OF POINTS IN THE PRISM
00198 !         IKLE HAS BEEN BUILT SO THAT IT COINCIDES WITH IKLE2D ON THE
00199 !         FIRST LAYER (SEE CPIKLE3).
00200 !
00201           II1=IKLE(IELEM,1)
00202           II2=IKLE(IELEM,2)
00203           II3=IKLE(IELEM,3)
00204           I1=II1+(IPLAN-1)*NPOIN2
00205           I2=II2+(IPLAN-1)*NPOIN2
00206           I3=II3+(IPLAN-1)*NPOIN2
00207           I4=I1+NPOIN2
00208           I5=I2+NPOIN2
00209           I6=I3+NPOIN2
00210 !
00211           X2=X(II2)-X(II1)
00212           X3=X(II3)-X(II1)
00213           Y2=Y(II2)-Y(II1)
00214           Y3=Y(II3)-Y(II1)
00215           SURFSUR3=XMUL*(X2*Y3-X3*Y2)/6.D0
00216 !
00217 !         EVERY TETRAHEDRON IN THE PRISM
00218 !
00219           DO I=1,3
00220 !
00221             IELEM3D=3*NELEM2*(IPLAN-1)+(I-1)*NELEM2+IELEM
00222 !
00223 !           CONTRIBUTION OF HORIZONTAL GRADIENTS
00224 !
00225             ALFA(1)=XMUL*H(IELEM3D,1)
00226             ALFA(2)=XMUL*H(IELEM3D,2)
00227             ALFA(3)=XMUL*H(IELEM3D,3)
00228             ALFA(4)=XMUL*H(IELEM3D,4)
00229 !
00230 !           CONTRIBUTION OF VERTICAL GRADIENTS
00231 !           ONLY FOR VERTICAL SEGMENTS
00232 !
00233 !           LOOP ON EVERY SEGMENT IN THE PRISM
00234 !
00235             DO ISEG=1,6
00236 !             SEE NUMBERING OF TETRAHEDRONS IN CPIKLE3
00237               IL1=ISEGT(ISEG,1)
00238               IL2=ISEGT(ISEG,2)
00239               J1=IKLE(IELEM3D,IL1)
00240               J2=IKLE(IELEM3D,IL2)
00241 !             VERTICAL SEGMENTS
00242               IF(J1.EQ.I1.AND.J2.EQ.I4) THEN
00243                 ALFA(IL1)=ALFA(IL1)-W(I1)*SURFSUR3
00244                 ALFA(IL2)=ALFA(IL2)+W(I1)*SURFSUR3
00245               ELSEIF(J1.EQ.I4.AND.J2.EQ.I1) THEN
00246                 ALFA(IL1)=ALFA(IL1)+W(I1)*SURFSUR3
00247                 ALFA(IL2)=ALFA(IL2)-W(I1)*SURFSUR3
00248               ELSEIF(J1.EQ.I2.AND.J2.EQ.I5) THEN
00249                 ALFA(IL1)=ALFA(IL1)-W(I2)*SURFSUR3
00250                 ALFA(IL2)=ALFA(IL2)+W(I2)*SURFSUR3
00251               ELSEIF(J1.EQ.I5.AND.J2.EQ.I2) THEN
00252                 ALFA(IL1)=ALFA(IL1)+W(I2)*SURFSUR3
00253                 ALFA(IL2)=ALFA(IL2)-W(I2)*SURFSUR3
00254               ELSEIF(J1.EQ.I3.AND.J2.EQ.I6) THEN
00255                 ALFA(IL1)=ALFA(IL1)-W(I3)*SURFSUR3
00256                 ALFA(IL2)=ALFA(IL2)+W(I3)*SURFSUR3
00257               ELSEIF(J1.EQ.I6.AND.J2.EQ.I3) THEN
00258                 ALFA(IL1)=ALFA(IL1)+W(I3)*SURFSUR3
00259                 ALFA(IL2)=ALFA(IL2)-W(I3)*SURFSUR3
00260               ENDIF
00261             ENDDO
00262 !
00263 !           N-SCHEME (SEE REFERENCE 2 PAGE 3, BOTTOM RIGHT)
00264 !
00265             SUMMAXK=MAX(0.D0,ALFA(1))
00266      &             +MAX(0.D0,ALFA(2))
00267      &             +MAX(0.D0,ALFA(3))
00268      &             +MAX(0.D0,ALFA(4))
00269             IF(SUMMAXK.GT.1.D-10) THEN
00270               SUMMAXK=-1.D0/SUMMAXK
00271               L12=SUMMAXK*MAX(0.D0,ALFA(1))*MIN(0.D0,ALFA(2))
00272               L13=SUMMAXK*MAX(0.D0,ALFA(1))*MIN(0.D0,ALFA(3))
00273               L14=SUMMAXK*MAX(0.D0,ALFA(1))*MIN(0.D0,ALFA(4))
00274               L23=SUMMAXK*MAX(0.D0,ALFA(2))*MIN(0.D0,ALFA(3))
00275               L24=SUMMAXK*MAX(0.D0,ALFA(2))*MIN(0.D0,ALFA(4))
00276               L34=SUMMAXK*MAX(0.D0,ALFA(3))*MIN(0.D0,ALFA(4))
00277               L21=SUMMAXK*MAX(0.D0,ALFA(2))*MIN(0.D0,ALFA(1))
00278               L31=SUMMAXK*MAX(0.D0,ALFA(3))*MIN(0.D0,ALFA(1))
00279               L41=SUMMAXK*MAX(0.D0,ALFA(4))*MIN(0.D0,ALFA(1))
00280               L32=SUMMAXK*MAX(0.D0,ALFA(3))*MIN(0.D0,ALFA(2))
00281               L42=SUMMAXK*MAX(0.D0,ALFA(4))*MIN(0.D0,ALFA(2))
00282               L43=SUMMAXK*MAX(0.D0,ALFA(4))*MIN(0.D0,ALFA(3))
00283             ELSE
00284 !             A SIMPLE SOLUTION, WITHOUT DIVISION
00285 !             PRINCIPLE : POINTS 1, 2, 3 GIVE THEIR CONTRIBUTION
00286 !             TO POINT 4, IT WORKS BECAUSE ALFA(4)=-ALFA(1)-ALFA(2)-ALFA(3)
00287 !             BUT IT IS NOT THE N-SCHEME
00288 !             BETTER THAN PUTTING 0 FOR MASS ERROR ?
00289 !             COULD BE USELESS DEPENDING ON FLUX CLIPPING AFTER
00290               L12 = 0.D0
00291               L13 = 0.D0
00292               L14 = MAX(ALFA(1),0.D0)
00293               L23 = 0.D0
00294               L24 = MAX(ALFA(2),0.D0)
00295               L34 = MAX(ALFA(3),0.D0)
00296               L21 = 0.D0
00297               L31 = 0.D0
00298               L41 = - MIN(ALFA(1),0.D0)
00299               L32 = 0.D0
00300               L42 = - MIN(ALFA(2),0.D0)
00301               L43 = - MIN(ALFA(3),0.D0)
00302             ENDIF
00303 !
00304             XM(01,IELEM3D) = L12
00305             XM(02,IELEM3D) = L13
00306             XM(03,IELEM3D) = L14
00307             XM(04,IELEM3D) = L23
00308             XM(05,IELEM3D) = L24
00309             XM(06,IELEM3D) = L34
00310             XM(07,IELEM3D) = L21
00311             XM(08,IELEM3D) = L31
00312             XM(09,IELEM3D) = L41
00313             XM(10,IELEM3D) = L32
00314             XM(11,IELEM3D) = L42
00315             XM(12,IELEM3D) = L43
00316 !
00317           ENDDO
00318 !
00319 !         END OF LOOP ON THE 3 TETRAHEDRONS
00320 !
00321         ENDDO
00322       ENDDO
00323 !
00324 !-----------------------------------------------------------------------
00325 !
00326 !     COMPUTES THE SUM OF EACH ROW (WITH A - SIGN)
00327 !     THE DIAGONAL TERMS ARE 0
00328 !
00329       IF(LEGO) THEN
00330 !
00331         DO IELEM = 1,NELEM
00332           T(IELEM,1) = -XM(01,IELEM)-XM(02,IELEM)-XM(03,IELEM)
00333           T(IELEM,2) = -XM(04,IELEM)-XM(05,IELEM)-XM(07,IELEM)
00334           T(IELEM,3) = -XM(06,IELEM)-XM(08,IELEM)-XM(10,IELEM)
00335           T(IELEM,4) = -XM(09,IELEM)-XM(11,IELEM)-XM(12,IELEM)
00336         ENDDO
00337 !
00338       ENDIF
00339 !
00340 !-----------------------------------------------------------------------
00341 !
00342       RETURN
00343       END

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