mt08pp.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE mt08pp
!                   *****************
!
     &( t,xm,xmul,sf,f,surfac,ikle,nelem,nelmax)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
!code
!+
!+                  /                     D
!+  A(I,J)=-XMUL   /  PSI2(J) *    F    * --( PSI1(I) ) D(OMEGA)
!+                /OMEGA                  DX
!+
!+  BEWARE THE MINUS SIGN !!
!+
!+  PSI1 AND PSI2: BASES OF TYPE P1 PRISM
!
!warning  THE JACOBIAN MUST BE POSITIVE
!warning  NEED TO CHECK THE SIGN; SEE USE IN DIFF3D!!!!!!!!!!!!!!!!
!
!history  J-M HERVOUET (LNH)    ; F  LEPEINTRE (LNH)
!+        28/11/94
!+        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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| F              |-->| FUNCTION USED IN THE FORMULA
!| FORMUL         |-->| FORMULA DESCRIBING THE RESULTING MATRIX
!| IKLE           |-->| CONNECTIVITY TABLE.
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| SF             |-->| STRUCTURE OF FUNCTIONS F
!| SURFAC         |-->| AREA OF 2D ELEMENTS
!| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
!| Z              |-->| ELEVATIONS OF POINTS
!| XM             |<->| OFF-DIAGONAL TERMS
!| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_mt08pp => mt08pp
!
      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
      DOUBLE PRECISION, INTENT(IN) :: F(*)
!
!     STRUCTURE OF F
!
      TYPE(bief_obj), INTENT(IN) :: SF
!
      DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
!
      DOUBLE PRECISION PZ1,XSU360
      DOUBLE PRECISION Q1,Q2,Q3,Q4,Q5,Q6
      DOUBLE PRECISION W14,W41,W25,W52,W63,W36
!
      INTEGER I1,I2,I3,I4,I5,I6,IELEM
!
!**********************************************************************
!
      xsu360 = xmul/360.d0
!
      IF(sf%ELM.NE.41) THEN
        WRITE(lu,1001) sf%ELM
1001    FORMAT(1x,'MT08PP (BIEF) : TYPE OF F NOT IMPLEMENTED: ',i6)
        CALL plante(1)
        stop
      ENDIF
!
!     LOOP ON THE ELEMENTS
!
      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)
!
        q1  =  f(i1)
        q2  =  f(i2)
        q3  =  f(i3)
        q4  =  f(i4)
        q5  =  f(i5)
        q6  =  f(i6)
!
!       INTERMEDIATE COMPUTATIONS
!
        pz1=-xsu360*surfac(ielem)
!
        w14 = q1+2*q4
        w41 = q4+2*q1
        w25 = q2+2*q5
        w52 = q5+2*q2
        w63 = q6+2*q3
        w36 = q3+2*q6
!
        t(ielem,1)=pz1*2*(3*w41+w52+w63)
        xm(ielem,18)=-t(ielem,1)
        xm(ielem,16)=pz1*(2*(w41+w52)+w63)
        xm(ielem,19)=-xm(ielem,16)
        xm(ielem,1) = xm(ielem,16)
        xm(ielem,22)=-xm(ielem,16)
        xm(ielem,2)=pz1*(2*(w41+w63)+w52)
        xm(ielem,20)=-xm(ielem,2)
        xm(ielem,17)= xm(ielem,2)
        xm(ielem,25)=-xm(ielem,2)
        t(ielem,2)=pz1*2*(w41+3*w52+w63)
        xm(ielem,23)= -t(ielem,2)
        xm(ielem,21)=pz1*(2*(w52+w63)+w41)
        xm(ielem,24)=-xm(ielem,21)
        xm(ielem,6) = xm(ielem,21)
        xm(ielem,26)=-xm(ielem,21)
        t(ielem,3)=pz1*2*(w41+w52+3*w63)
        xm(ielem,27)=-t(ielem,3)
        xm(ielem,3)=pz1*2*(3*w14+w25+w36)
        t(ielem,4)=-xm(ielem,3)
        xm(ielem,7)=pz1*(2*(w14+w25)+w36)
        xm(ielem,28)=-xm(ielem,7)
        xm(ielem,4) = xm(ielem,7)
        xm(ielem,13)=-xm(ielem,7)
        xm(ielem,10)=pz1*(2*(w14+w36)+w25)
        xm(ielem,29)=-xm(ielem,10)
        xm(ielem,5) = xm(ielem,10)
        xm(ielem,14)=-xm(ielem,10)
        xm(ielem,8)=pz1*2*(w14+3*w25+w36)
        t(ielem,5)=-xm(ielem,8)
        xm(ielem,11)=pz1*(2*(w25+w36)+w14)
        xm(ielem,30)=-xm(ielem,11)
        xm(ielem,9) = xm(ielem,11)
        xm(ielem,15)=-xm(ielem,11)
        xm(ielem,12)=pz1*2*(w14+w25+3*w36)
        t(ielem,6)=-xm(ielem,12)
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END