mt08tt.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE mt08tt
!                   *****************
!
     &( t,xm,xmul,x,y,sf,f,ikle,nelem,nelmax)
!
!***********************************************************************
! BIEF   V6P1                                  21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
!code
!+
!+               /
!+  A    = XMUL / P   F . GRAD(P ) * J(X,Y) DXDY
!+   I J       /S  J            I
!+
!+  BY ELEMENTARY CELL; THE ELEMENT IS THE P1 TRIANGLE
!+
!+  J(X,Y): JACOBIAN OF THE ISOPARAMETRIC TRANSFORMATION
!
!note     ONLY THE Z COMPONENT IS TREATED HERE !!
!
!warning  THE JACOBIAN MUST BE POSITIVE
!
!history  J-M HERVOUET (LNH)
!+        21/03/02
!+        V5P3
!+
!
!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
!| XM             |<->| OFF-DIAGONAL TERMS
!| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_mt08tt => mt08tt
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,INTENT(IN) :: NELEM,NELMAX
      INTEGER,INTENT(IN) :: IKLE(nelmax,4)
!
      DOUBLE PRECISION,INTENT(INOUT) :: T(nelmax,4),XM(nelmax,12)
!
      DOUBLE PRECISION,INTENT(IN) :: XMUL
      DOUBLE PRECISION,INTENT(IN) :: F(*),X(*),Y(*)
!
!     STRUCTURE OF F
!
      TYPE(bief_obj),INTENT(IN) :: SF
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
!
      DOUBLE PRECISION X2,X3,X4,Y2,Y3,Y4,F1,F2,F3,F4,XSUR120
!
      INTEGER I1,I2,I3,I4,IELEM
!
!**********************************************************************
!
      xsur120 = xmul/120.d0
!
      IF(sf%ELM.NE.31.AND.sf%ELM.NE.51) THEN
        WRITE(lu,1001) sf%ELM
1001    FORMAT(1x,'MT08TT (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)
!
        x2 = x(i2)-x(i1)
        x3 = x(i3)-x(i1)
        x4 = x(i4)-x(i1)
!
        y2 = y(i2)-y(i1)
        y3 = y(i3)-y(i1)
        y4 = y(i4)-y(i1)
!
        f1 = f(i1)
        f2 = f(i2)
        f3 = f(i3)
        f4 = f(i4)
!
        t(ielem,1)=(
     & x2*y4*f3-x3*y4*f3+x4*y3*f2-y2*x4*f2+x2*y4*f2+2*y2*x3*f1
     &-2*x3*y4*f1+y2*x3*f2+y2*x3*f3+2*x4*y3*f1+2*x2*y4*f1-2*y2*x4*f1
     &-2*x2*y3*f1-y2*x4*f3-x2*y3*f2-x2*y3*f3+x4*y3*f3-x3*y4*f4
     &+x2*y4*f4-y2*x4*f4-x2*y3*f4+y2*x3*f4+x4*y3*f4-x3*y4*f2 )*xsur120
!
        t(ielem,2)=(
     & x3*y4*f3-2*x4*y3*f2+x3*y4*f1-x4*y3*f1-x4*y3*f3+x3*y4*f4
     &-x4*y3*f4+2*x3*y4*f2                        )*xsur120
!
        t(ielem,3)=(-x2*y4+y2*x4)*(f1+2*f3+f2+f4) *xsur120
!
        t(ielem,4)=(
     & -y2*x3*f1-y2*x3*f2-y2*x3*f3+x2*y3*f1+x2*y3*f2+x2*y3*f3+
     &2*x2*y3*f4-2*y2*x3*f4                       )*xsur120
!
        xm(ielem,01)=(
     & x2*y4*f3-x3*y4*f3+2*x4*y3*f2-2*y2*x4*f2+2*x2*y4*f2+y2*x3*f1
     &-x3*y4*f1+2*y2*x3*f2+y2*x3*f3+x4*y3*f1+x2*y4*f1-y2*x4*f1-x2*y3
     &*f1-y2*x4*f3-2*x2*y3*f2-x2*y3*f3+x4*y3*f3
     &-x3*y4*f4+x2*y4*f4-y2*x4*f4
     &-x2*y3*f4+y2*x3*f4+x4*y3*f4-2*x3*y4*f2    )*xsur120
!
        xm(ielem,02)=(
     &-x3*y4+x4*y3+x2*y4-y2*x4-x2*y3+y2*x3)*(f1+2*f3+f2+f4)*xsur120
!
        xm(ielem,03)=(
     & x2*y4*f3-x3*y4*f3+x4*y3*f2-y2*x4*f2+x2*y4*f2+y2*x3*f1-x3*y4*f1
     &+y2*x3*f2+y2*x3*f3+x4*y3*f1+x2*y4*f1-y2*x4*f1-x2*y3*f1-y2*x4*f3
     &-x2*y3*f2-x2*y3*f3+x4*y3*f3-2*x3*y4*f4+2*x2*y4*f4-2*y2*x4*f4-2
     &*x2*y3*f4+2*y2*x3*f4+2*x4*y3*f4-x3*y4*f2)*xsur120
!
        xm(ielem,04)= -(-x3*y4+x4*y3)*(f1+2*f3+f2+f4)*xsur120
!
        xm(ielem,05)=( x3*y4*f3-x4*y3*f2+x3*y4*f1-x4*y3*f1-x4*y3*f3
     &                 +2*x3*y4*f4-2*x4*y3*f4+x3*y4*f2)*xsur120
!
        xm(ielem,06)=( -x2*y4*f3+y2*x4*f2-x2*y4*f2-x2*y4*f1+y2*x4*f1
     &                 +y2*x4*f3-2*x2*y4*f4+2*y2*x4*f4)*xsur120
!
        xm(ielem,07)=( x3*y4*f3-x4*y3*f2+2*x3*y4*f1-2*x4*y3*f1
     &                -x4*y3*f3+x3*y4*f4-x4*y3*f4+x3*y4*f2)*xsur120
!
        xm(ielem,08)=( -x2*y4*f3+y2*x4*f2-x2*y4*f2-2*x2*y4*f1
     &                 +2*y2*x4*f1+y2*x4*f3-x2*y4*f4+y2*x4*f4)*xsur120
!
        xm(ielem,09)=( -x2*y4*f3+2*y2*x4*f2-2*x2*y4*f2-x2*y4*f1
     &                 +y2*x4*f1+y2*x4*f3-x2*y4*f4+y2*x4*f4 )*xsur120
!
        xm(ielem,10)=( -2*y2*x3*f1-y2*x3*f2-y2*x3*f3+2*x2*y3*f1
     &                 +x2*y3*f2+x2*y3*f3+x2*y3*f4-y2*x3*f4)*xsur120
!
        xm(ielem,11)=(-y2*x3*f1-2*y2*x3*f2-y2*x3*f3+x2*y3*f1
     &                +2*x2*y3*f2+x2*y3*f3+x2*y3*f4-y2*x3*f4)*xsur120
!
        xm(ielem,12)= -(-x2*y3+y2*x3)*(f1+2*f3+f2+f4)*xsur120
!
      ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
      RETURN
      END