gettriebe.f

Go to the documentation of this file.

!                   ********************
                    SUBROUTINE gettriebe
!                   ********************
!
     &(xaux,ad,ax,teta,ikle,npoin,nelem,nelmax,mesh,ielm3,nelem2,nplan,
     & knolg)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    GETS THE TRIDIAGONAL PART OF A DIFFUSION MATRIX ON
!+                 PRISMS AND REMOVES IT FROM THE INITIAL MATRIX.
!code
!+            IF MTRI IS THIS TRIDIAGONAL PART, MAUX THE RESULT AND MDIF
!+            THE DIFFUSION MATRIX, THIS SUBROUTINE DOES:
!+
!+            MAUX = TETA * MTRI
!+            MDIF CHANGED INTO (1-TETA) * MDIF
!
!warning  THE JACOBIAN MUST BE POSITIVE
!
!history  J-M HERVOUET (LNHE)
!+        13/08/08
!+        V5P9
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (LNHE)
!+        29/11/2011
!+        V6P2
!+   Element 51 programmed, KNOLG added
!+
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AD             |-->| DIAGONAL TERMS OF MATRIX
!| AX             |-->| OFF-DIAGONAL TERMS OF MATRIX
!| IELM3          |-->| TYPE OF ELEMENT
!| IKLE           |-->| CONNECTIVITY TABLE
!| KNOLG          |-->| ORIGINAL (I.E. SCALAR MODE) GLOBAL NUMBER OF POINTS
!| MESH           |-->| MESH STRUCTURE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELEM2         |-->| NUMBER OF TRIANGLES OF ORIGINAL 2D MESH
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPLAN          |-->| NUMBER OF PLANES IN THE ORIGINAL MESH
!| NPOIN          |-->| NUMBER OF POINTS
!| TETA           |-->| COEFFICIENT USED IN THE RESULT
!| XAUX           |<--| THE RESULTING MATRIX
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_gettriebe => gettriebe
      USE declarations_telemac, ONLY : tetra,isegt
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN,IELM3,NELEM2,NPLAN
      INTEGER, INTENT(IN) :: IKLE(nelmax,*),KNOLG(npoin)
!
      DOUBLE PRECISION, INTENT(IN)    :: TETA
      DOUBLE PRECISION, INTENT(INOUT) :: XAUX(npoin,*),AX(nelmax,*)
      DOUBLE PRECISION, INTENT(INOUT) :: AD(npoin)
!
      TYPE(bief_mesh) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I1,I2,I3,I4,I5,I6,IELEM,IAN,ICOM,NPOIN2,IPLAN,K
      INTEGER S1,S2,S3,IT1,IT2,IELEM3D,L1,L2,ISEG
!
!     TETRA : WILL GIVE THE LOCAL NUMBERS OF POINTS IN THE PRISM
!             THE 0 CORRESPOND TO SITUATIONS
!             THAT NEVER HAPPEN (TETRA(1,1,1,... OR TETRA(2,2,2,...)
!             SEE ALSO CPIKLE3 AND FLUX_EF_VF_3D, WITH SIMILAR USE
!     INTEGER TETRA(2,2,2,3,4)
!     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,
!    &             0,2,2,2,2,2,2,0,0,6,6,6,6,6,6,0,0,3,1,2,2,1,3,0,
!    &             0,3,3,3,3,3,3,0,0,5,5,5,5,5,5,0,0,2,3,4,1,6,5,0,
!    &             0,4,5,4,6,6,5,0,0,2,3,3,1,2,1,0,0,4,5,3,6,2,1,0 /
!     EBE SYMMETRIC STORAGE OF OFF-DIAGONAL TERMS FOR TETRAHEDRONS
!     0 ARE PUT FOR DIAGONALS
      INTEGER :: STO(4,4)
      parameter( sto = reshape( (/
     &           0, 1, 2, 3,
     &           1, 0, 4, 5,
     &           2, 4, 0, 6,
     &           3, 5, 6, 0 /), shape=(/ 4,4 /) ) )
!     CORRESPONDING OFF-DIAGONAL TERMS
!     DATA STO /1-1,2-1,3-1,4-1,
!    &          1-2,2-2,3-2,4-2,
!    &          1-3,2-3,3-3,4-3,
!    &          1-4,2-4,3-4,4-4/
!     INTEGER ISEGT(6,2)
!     DATA ISEGT/1,2,3,1,2,3,2,3,1,4,4,4/
!
!-----------------------------------------------------------------------
!
!     CONSIDERS HERE THAT NPOIN < NELMAX TO USE XAUX AS XAUX(NPOIN,3)
!
!     XAUX(I,1) IS COEFFICIENT OF POINT BELOW I IN EQUATION OF POINT I
!     XAUX(I,2) IS THE DIAGONAL
!     XAUX(I,3) IS COEFFICIENT OF POINT ABOVE I IN EQUATION OF POINT I
!
!-----------------------------------------------------------------------
!     INITIALISES THE DIAGONAL AND OFF-DIAGONAL TERMS
!-----------------------------------------------------------------------
!
      CALL ov('X=C     ', x=xaux(1,1), c=0.d0, dim1=npoin)
      CALL ov('X=CY    ', x=xaux(1,2), y=ad, c=teta, dim1=npoin)
      CALL ov('X=C     ', x=xaux(1,3), c=0.d0, dim1=npoin)
!
      CALL ov('X=CX    ',x=ad, c=1.d0-teta, dim1=npoin)
!
!-----------------------------------------------------------------------
!     ADDS TRIDIAGONAL TERMS
!-----------------------------------------------------------------------
!
      IF(ielm3.EQ.41) THEN
!
        DO ielem=1,nelem
!
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
          i4=ikle(ielem,4)
          i5=ikle(ielem,5)
          i6=ikle(ielem,6)
          xaux(i1,3)=xaux(i1,3)+teta*ax(ielem,03) ! TERM 1-4
          xaux(i2,3)=xaux(i2,3)+teta*ax(ielem,08) ! TERM 2-5
          xaux(i3,3)=xaux(i3,3)+teta*ax(ielem,12) ! TERM 3-6
          xaux(i4,1)=xaux(i4,1)+teta*ax(ielem,03) ! TERM 4-1
          xaux(i5,1)=xaux(i5,1)+teta*ax(ielem,08) ! TERM 5-2
          xaux(i6,1)=xaux(i6,1)+teta*ax(ielem,12) ! TERM 6-3
!
          ax(ielem,03)=ax(ielem,03)*(1.d0-teta)
          ax(ielem,08)=ax(ielem,08)*(1.d0-teta)
          ax(ielem,12)=ax(ielem,12)*(1.d0-teta)
!
        ENDDO
!
      ELSEIF(ielm3.EQ.51) THEN
!
        DO iplan=1,nplan-1
          DO ielem=1,nelem2
!
!           HERE LOWER LEVEL OF ELEMENTS ALWAYS TAKEN
!           TO FIND THE WAY THE PRISM HAS BEEN CUT BY LOOKING
!           AT GLOBAL NUMBERS OF POINTS
!           THIS PART IS THUS COMMON TO ALL PLANES
!           IKLE 3D IS TAKEN, COULD BE IKLE 2D AS WELL
            i1=ikle(ielem,1)
            i2=ikle(ielem,2)
            i3=ikle(ielem,3)
            IF(ncsize.GT.1) THEN
              i1=knolg(i1)
              i2=knolg(i2)
              i3=knolg(i3)
            ENDIF
!           THIS IS DONE LIKE IN CPIKLE3 TO USE ARRAY TETRA
            IF(i1.GT.i2) THEN
              s1=1
            ELSE
              s1=2
            ENDIF
            IF(i2.GT.i3) THEN
              s2=1
            ELSE
              s2=2
            ENDIF
            IF(i3.GT.i1) THEN
              s3=1
            ELSE
              s3=2
            ENDIF
!
!           NOW TAKING CONTRIBUTIONS OF TETRAHEDRON K= 1, 2 AND 3
!
            DO k=1,3
              ielem3d=3*(iplan-1)*nelem2+ielem+(k-1)*nelem2
!             SEGMENTS 1 TO 6
              DO iseg=1,6
!               LOCAL NUMBERS OF 2 POINTS OF SEGMENT IN THE TETRAHEDRON
                l1=isegt(iseg,1)
                l2=isegt(iseg,2)
!               GLOBAL NUMBERS OF 2 POINTS OF SEGMENT
                i1=ikle(ielem3d,l1)
                i2=ikle(ielem3d,l2)
!               NUMBERS OF 2 POINTS OF SEGMENT IN THE ORIGINAL PRISM
                it1=tetra(s1,s2,s3,k,l1)
                it2=tetra(s1,s2,s3,k,l2)
!               WE LOOK FOR VERTICALS OF THE ORIGINAL PRISM
                IF(it1.EQ.1.AND.it2.EQ.4) THEN
                  xaux(i1,3)=xaux(i1,3)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 1-4
                  xaux(i2,1)=xaux(i2,1)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 4-1
                  ax(ielem3d,sto(l1,l2))=
     &            ax(ielem3d,sto(l1,l2))*(1.d0-teta)
                ELSEIF(it1.EQ.4.AND.it2.EQ.1) THEN
                  xaux(i1,1)=xaux(i1,1)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 4-1
                  xaux(i2,3)=xaux(i2,3)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 1-4
                  ax(ielem3d,sto(l1,l2))=
     &            ax(ielem3d,sto(l1,l2))*(1.d0-teta)
                ELSEIF(it1.EQ.2.AND.it2.EQ.5) THEN
                  xaux(i1,3)=xaux(i1,3)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 2-5
                  xaux(i2,1)=xaux(i2,1)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 5-2
                  ax(ielem3d,sto(l1,l2))=
     &            ax(ielem3d,sto(l1,l2))*(1.d0-teta)
                ELSEIF(it1.EQ.5.AND.it2.EQ.2) THEN
                  xaux(i1,1)=xaux(i1,1)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 5-2
                  xaux(i2,3)=xaux(i2,3)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 2-5
                  ax(ielem3d,sto(l1,l2))=
     &            ax(ielem3d,sto(l1,l2))*(1.d0-teta)
                ELSEIF(it1.EQ.3.AND.it2.EQ.6) THEN
                  xaux(i1,3)=xaux(i1,3)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 3-6
                  xaux(i2,1)=xaux(i2,1)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 6-3
                  ax(ielem3d,sto(l1,l2))=
     &            ax(ielem3d,sto(l1,l2))*(1.d0-teta)
                ELSEIF(it1.EQ.6.AND.it2.EQ.3) THEN
                  xaux(i1,1)=xaux(i1,1)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 6-3
                  xaux(i2,3)=xaux(i2,3)+teta*ax(ielem3d,sto(l1,l2)) ! TERM 3-6
                  ax(ielem3d,sto(l1,l2))=
     &            ax(ielem3d,sto(l1,l2))*(1.d0-teta)
                ENDIF
              ENDDO
            ENDDO
!
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
      ELSE
        WRITE(lu,*) 'GETTRIEBE: ELEMENT NOT IMPLEMENTED:',ielm3
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     PARALLEL MODE
!
      IF(ncsize.GT.1) THEN
        ian    = 3
        icom   = 2
        npoin2 = bief_nbpts(11,mesh)
        CALL parcom2(xaux(1,1),xaux(1,2),xaux(1,3),
     &               npoin2,nplan,icom,ian,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END