mt14tt.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE mt14tt
!                   *****************
!
     &(t,xm,lego,xmul,sw,w,h,
     & x,y,ikle,nelem,nelmax,nplan,npoin2)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    BUILDS COEFFICIENTS LAMBDA(I,J) FOR N-TYPE MURD SCHEME.
!+
!+            THE ELEMENT IS THE P1 TETRAHEDRON.
!
!warning  THE JACOBIAN MUST BE POSITIVE
!+        Only for prisms cut into tetrahedra so far.
!
!reference  J-M HERVOUET THESIS OR BOOK: MURD SCHEME IN 3 DIMENSIONS
!+             CHAPTER 6.
!+
!+          Computational Fluid dynamics '92 volume 1. Multidimensional
!+          upwind schemes for scalar advection on tetrahedral meshes.
!+          G. Bourgois, H. Deconinck, P.L. Roe and R. Struijs
!
!history  J-M HERVOUET (LNHE)
!+        18/09/2011
!+        V6P1
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| H              |-->| FUNCTION USED IN THE FORMULA
!| IKLE           |-->| CONNECTIVITY TABLE
!| LEGO           |-->| LOGICAL. IF YES RESULT ASSEMBLED.
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELEM2         |-->| NUMBER OF TRIANGLES IN THE UNDERLYING 2D MESH
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPLAN          |-->| NUMBER OF PLANES IN THE MESH
!| NPOIN2         |-->| NUMBER OF POINTS IN THE UNDERLYING 2D MESH
!| SW             |-->| BIEF_OBJ STRUCTURE OF W
!| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
!| W              |-->| FUNCTION USED IN THE FORMULA
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH
!| Y              |-->| ORDINATES OF POINTS IN THE MESH
!| XM             |<->| OFF-DIAGONAL TERMS
!| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief           !, EX_MT14TT => MT14TT
      USE declarations_telemac, ONLY : isegt
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,NPLAN,NPOIN2
      INTEGER, INTENT(IN) :: IKLE(nelmax,4)
!
      DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,4),XM(12,nelmax)
!
      DOUBLE PRECISION, INTENT(IN) :: XMUL
      DOUBLE PRECISION, INTENT(IN) :: W(*)
      DOUBLE PRECISION, INTENT(IN) :: H(nelmax,4)
!
      LOGICAL, INTENT(IN) :: LEGO
!
!     STRUCTURES DE U,V,W
!
      TYPE(bief_obj), INTENT(IN) :: SW
!
      DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
!
      INTEGER IELEM,IPLAN,J1,J2,ISEG,II1,II2,II3,IELEM3D,NELEM2
!
      DOUBLE PRECISION X2,X3,Y2,Y3
      DOUBLE PRECISION L12,L13,L14,L23,L24,L34,L21,L31,L41,L32,L42,L43
      DOUBLE PRECISION ALFA(4),SURFSUR3,SUMMAXK
!
      INTEGER IL1,IL2,I1,I2,I3,I4,I5,I6,I
!
!-----------------------------------------------------------------------
!
!     THE SIX SEGMENTS IN A TETRAHEDRON
!     ISEGT(ISEG,1 OR 2) : FIRST OR SECOND POINT OF SEGMENT ISEG
!     INTEGER ISEGT(6,2)
!     DATA ISEGT/1,2,3,1,2,3,2,3,1,4,4,4/
!
!=======================================================================
!
!     RETRIEVING THE NUMBER OF ELEMENTS THE 2D MESH
!
      nelem2=nelmax/3/(nplan-1)
!                                                         *   *
!     COMPUTES THE COEFFICIENTS A(I) (INTEGRAL OF H U.GRAD(PSI), ONLY
!     CONSIDERING THE HORIZONTAL PART OF U). SEE JMH THESIS
!
!     NOTE: THIS VC04TT CALL IS ALREADY DONE BY FLUX3D IN TELEMAC-3D
!           THROUGH A CALL TO VECTOR. THE EQUIVALENT OF T HAS BEEN
!           SAVED IN WHAT IS HERE H
!
!     FORMUL='             HOR'
!     CALL VC04TT(XMUL,SU,SV,SW,U,V,W,SF,SG,SH,F,G,H,X,Y,Z,
!    * IKLE(1,1),IKLE(1,2),IKLE(1,3),IKLE(1,4),
!    * NELEM,NELMAX,T(1,1),T(1,2),T(1,3),T(1,4),
!    * SPECAD,FORMUL,NPLAN)
!
!     NOW T IS RETRIEVED BY STORAGE DONE INTO FLUX3D
!     FUNCTION H MUST CONTAIN THE EQUIVALENT OF T ABOVE, WHICH IS
!     THE NON ASSEMBLED VALUE OF FLUINT
!
!     DO IELEM=1,NELEM
!       T(IELEM,1)=H(         IELEM)
!       T(IELEM,2)=H(  NELMAX+IELEM)
!       T(IELEM,3)=H(2*NELMAX+IELEM)
!       T(IELEM,4)=H(3*NELMAX+IELEM)
!     ENDDO
!
!     FORMUL NORMALLY STARTS WITH VGRADP OR VGRADP2, OR VGRADP22 TO KNOW
!     SIGMAG AND SPECAD, USELESS HERE.
!
!-----------------------------------------------------------------------
!
      IF(sw%ELM.NE.51) THEN
        WRITE(lu,*) 'MT14TT DISCRETISATION NON PREVUE'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      DO iplan=1,nplan-1
        DO ielem=1,nelem2
!
!         THE SIX GLOBAL NUMBERS OF POINTS IN THE PRISM
!         IKLE HAS BEEN BUILT SO THAT IT COINCIDES WITH IKLE2D ON THE
!         FIRST LAYER (SEE CPIKLE3).
!
          ii1=ikle(ielem,1)
          ii2=ikle(ielem,2)
          ii3=ikle(ielem,3)
          i1=ii1+(iplan-1)*npoin2
          i2=ii2+(iplan-1)*npoin2
          i3=ii3+(iplan-1)*npoin2
          i4=i1+npoin2
          i5=i2+npoin2
          i6=i3+npoin2
!
          x2=x(ii2)-x(ii1)
          x3=x(ii3)-x(ii1)
          y2=y(ii2)-y(ii1)
          y3=y(ii3)-y(ii1)
          surfsur3=xmul*(x2*y3-x3*y2)/6.d0
!
!         EVERY TETRAHEDRON IN THE PRISM
!
          DO i=1,3
!
            ielem3d=3*nelem2*(iplan-1)+(i-1)*nelem2+ielem
!
!           CONTRIBUTION OF HORIZONTAL GRADIENTS
!
            alfa(1)=xmul*h(ielem3d,1)
            alfa(2)=xmul*h(ielem3d,2)
            alfa(3)=xmul*h(ielem3d,3)
            alfa(4)=xmul*h(ielem3d,4)
!
!           CONTRIBUTION OF VERTICAL GRADIENTS
!           ONLY FOR VERTICAL SEGMENTS
!
!           LOOP ON EVERY SEGMENT IN THE PRISM
!
            DO iseg=1,6
!             SEE NUMBERING OF TETRAHEDRONS IN CPIKLE3
              il1=isegt(iseg,1)
              il2=isegt(iseg,2)
              j1=ikle(ielem3d,il1)
              j2=ikle(ielem3d,il2)
!             VERTICAL SEGMENTS
              IF(j1.EQ.i1.AND.j2.EQ.i4) THEN
                alfa(il1)=alfa(il1)-w(i1)*surfsur3
                alfa(il2)=alfa(il2)+w(i1)*surfsur3
              ELSEIF(j1.EQ.i4.AND.j2.EQ.i1) THEN
                alfa(il1)=alfa(il1)+w(i1)*surfsur3
                alfa(il2)=alfa(il2)-w(i1)*surfsur3
              ELSEIF(j1.EQ.i2.AND.j2.EQ.i5) THEN
                alfa(il1)=alfa(il1)-w(i2)*surfsur3
                alfa(il2)=alfa(il2)+w(i2)*surfsur3
              ELSEIF(j1.EQ.i5.AND.j2.EQ.i2) THEN
                alfa(il1)=alfa(il1)+w(i2)*surfsur3
                alfa(il2)=alfa(il2)-w(i2)*surfsur3
              ELSEIF(j1.EQ.i3.AND.j2.EQ.i6) THEN
                alfa(il1)=alfa(il1)-w(i3)*surfsur3
                alfa(il2)=alfa(il2)+w(i3)*surfsur3
              ELSEIF(j1.EQ.i6.AND.j2.EQ.i3) THEN
                alfa(il1)=alfa(il1)+w(i3)*surfsur3
                alfa(il2)=alfa(il2)-w(i3)*surfsur3
              ENDIF
            ENDDO
!
!           N-SCHEME (SEE REFERENCE 2 PAGE 3, BOTTOM RIGHT)
!
            summaxk=max(0.d0,alfa(1))
     &             +max(0.d0,alfa(2))
     &             +max(0.d0,alfa(3))
     &             +max(0.d0,alfa(4))
            IF(summaxk.GT.1.d-10) THEN
              summaxk=-1.d0/summaxk
              l12=summaxk*max(0.d0,alfa(1))*min(0.d0,alfa(2))
              l13=summaxk*max(0.d0,alfa(1))*min(0.d0,alfa(3))
              l14=summaxk*max(0.d0,alfa(1))*min(0.d0,alfa(4))
              l23=summaxk*max(0.d0,alfa(2))*min(0.d0,alfa(3))
              l24=summaxk*max(0.d0,alfa(2))*min(0.d0,alfa(4))
              l34=summaxk*max(0.d0,alfa(3))*min(0.d0,alfa(4))
              l21=summaxk*max(0.d0,alfa(2))*min(0.d0,alfa(1))
              l31=summaxk*max(0.d0,alfa(3))*min(0.d0,alfa(1))
              l41=summaxk*max(0.d0,alfa(4))*min(0.d0,alfa(1))
              l32=summaxk*max(0.d0,alfa(3))*min(0.d0,alfa(2))
              l42=summaxk*max(0.d0,alfa(4))*min(0.d0,alfa(2))
              l43=summaxk*max(0.d0,alfa(4))*min(0.d0,alfa(3))
            ELSE
!             A SIMPLE SOLUTION, WITHOUT DIVISION
!             PRINCIPLE : POINTS 1, 2, 3 GIVE THEIR CONTRIBUTION
!             TO POINT 4, IT WORKS BECAUSE ALFA(4)=-ALFA(1)-ALFA(2)-ALFA(3)
!             BUT IT IS NOT THE N-SCHEME
!             BETTER THAN PUTTING 0 FOR MASS ERROR ?
!             COULD BE USELESS DEPENDING ON FLUX CLIPPING AFTER
              l12 = 0.d0
              l13 = 0.d0
              l14 = max(alfa(1),0.d0)
              l23 = 0.d0
              l24 = max(alfa(2),0.d0)
              l34 = max(alfa(3),0.d0)
              l21 = 0.d0
              l31 = 0.d0
              l41 = - min(alfa(1),0.d0)
              l32 = 0.d0
              l42 = - min(alfa(2),0.d0)
              l43 = - min(alfa(3),0.d0)
            ENDIF
!
            xm(01,ielem3d) = l12
            xm(02,ielem3d) = l13
            xm(03,ielem3d) = l14
            xm(04,ielem3d) = l23
            xm(05,ielem3d) = l24
            xm(06,ielem3d) = l34
            xm(07,ielem3d) = l21
            xm(08,ielem3d) = l31
            xm(09,ielem3d) = l41
            xm(10,ielem3d) = l32
            xm(11,ielem3d) = l42
            xm(12,ielem3d) = l43
!
          ENDDO
!
!         END OF LOOP ON THE 3 TETRAHEDRONS
!
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
!     COMPUTES THE SUM OF EACH ROW (WITH A - SIGN)
!     THE DIAGONAL TERMS ARE 0
!
      IF(lego) THEN
!
        DO ielem = 1,nelem
          t(ielem,1) = -xm(01,ielem)-xm(02,ielem)-xm(03,ielem)
          t(ielem,2) = -xm(04,ielem)-xm(05,ielem)-xm(07,ielem)
          t(ielem,3) = -xm(06,ielem)-xm(08,ielem)-xm(10,ielem)
          t(ielem,4) = -xm(09,ielem)-xm(11,ielem)-xm(12,ielem)
        ENDDO
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END