mt01pp.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE mt01pp
!                   *****************
!
     &( t,xm,xmul,z,surfac,ikle,nelem,nelmax)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    BUILDS THE MASS MATRIX FOR P1 PRISMS.
!code
!+    COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
!+
!+                                 /
!+                    A    = XMUL /  (P *P )*J(X,Y) DXDY
!+                     I J       /S    I  J
!+
!+    BY ELEMENTARY CELL (REAL MESH)
!+
!+    J(X,Y) : JACOBIAN OF THE ISOPARAMETRIC TRANSFORMATION
!
!warning  THE JACOBIAN MUST BE POSITIVE
!
!history  J-M HERVOUET (LNH)    ; F  LEPEINTRE (LNH)
!+        21/05/2010
!+        V6P0
!+
!
!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  U.H.MErkel
!+        18/07/2012
!+        V6P2
!+   Replaced EPSILON with CHOUIA due to nag compiler problems
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IKLE           |-->| CONNECTIVITY TABLE.
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| SURFAC         |-->| AREA OF TRIANGLES
!| T              |<--| NON ASSEMBLED DIAGONAL
!| XM             |<--| NON ASSEMBLED OFF-DIAGONAL TERMS
!| XMUL           |-->| MULTIPLICATION FACTOR
!| Z              |-->| ELEVATIONS OF POINTS IN THE MESH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_mt01pp => mt01pp
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NELEM,NELMAX
      INTEGER, INTENT(IN)             :: IKLE(nelmax,6)
      DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(nelmax,30)
      DOUBLE PRECISION, INTENT(IN)    :: XMUL,Z(*),SURFAC(nelmax)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
!
      INTEGER IELEM
      DOUBLE PRECISION SUR360,COEF,H1,H2,H3,HT
!
      DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-3
!
!-----------------------------------------------------------------------
!
      sur360 = xmul / 360.d0
!
!   LOOP ON THE ELEMENTS
!
      DO ielem = 1,nelem
!
        coef = surfac(ielem) * sur360
!
!       TREATS HERE THE DRY ZONES
!
        h1 = max(z(ikle(ielem,4)) - z(ikle(ielem,1)),chouia) * coef
        h2 = max(z(ikle(ielem,5)) - z(ikle(ielem,2)),chouia) * coef
        h3 = max(z(ikle(ielem,6)) - z(ikle(ielem,3)),chouia) * coef
        ht = h1 + h2 + h3
!
!-----------------------------------------------------------------------
!
!  EXTRA-DIAGONAL TERMS
!
        xm(ielem,4)  = h1 + h2 + ht
        xm(ielem,5)  = h1 + h3 + ht
        xm(ielem,9)  = h2 + h3 + ht
        xm(ielem,7)  = xm(ielem,4)
        xm(ielem,10) = xm(ielem,5)
        xm(ielem,11) = xm(ielem,9)
!
        xm(ielem,3)  =  4*h1 + ht + ht
        xm(ielem,8)  =  4*h2 + ht + ht
        xm(ielem,12) =  4*h3 + ht + ht
!
        xm(ielem,1)  = xm(ielem,4) + xm(ielem,4)
        xm(ielem,2)  = xm(ielem,5) + xm(ielem,5)
        xm(ielem,6)  = xm(ielem,9) + xm(ielem,9)
        xm(ielem,13) = xm(ielem,1)
        xm(ielem,14) = xm(ielem,2)
        xm(ielem,15) = xm(ielem,6)
!
!  DIAGONAL TERMS
!
        t(ielem,1) = xm(ielem,3)  + xm(ielem,3)
        t(ielem,2) = xm(ielem,8)  + xm(ielem,8)
        t(ielem,3) = xm(ielem,12) + xm(ielem,12)
        t(ielem,4) = t(ielem,1)
        t(ielem,5) = t(ielem,2)
        t(ielem,6) = t(ielem,3)
!
!   END OF THE LOOP ON THE ELEMENTS
!
      ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
      RETURN
      END