mt12aa.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE mt12aa
!                   *****************
!
     &(  a11 , a12 , a13 ,
     &   a21 , a22 , a23 ,
     &   a31 , a32 , a33 ,
     &   xmul,sf,su,sv,f,u,v,
     &   xel,yel,surdet,
     &   ikle1,ikle2,ikle3,
     &   nelem,nelmax,icoord)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
!code
!+  EXAMPLE WITH ICOORD=1
!+
!+                  /                -->  --->
!+  A(I,J) = XMUL  /   PSIJ GRAD(F)   U .GRAD(PSII) D(OMEGA)
!+                /OMEGA
!+                                 -->
!+                      F VECTOR    U  VECTOR WITH COMPONENTS U,V,W
!+
!+
!+  BEWARE THE MINUS SIGN !!
!+
!+  PSI1: LINEAR
!+  PSI2: LINEAR
!+
!+  IT WOULD BE A DERIVATIVE WRT Y WITH ICOORD=2
!+  IT WOULD BE A DERIVATIVE WRT Z WITH ICOORD=3
!
!warning  THE JACOBIAN MUST BE POSITIVE
!
!history  J-M HERVOUET (LNH)    ; C MOULIN (LNH)
!+        09/12/94
!+        V5P1
!+
!
!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)
!+        15/06/2012
!+        V6P2
!+   SURFAC changed into SURDET and formulas changed accordingly
!+   This is just optimisation.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A11            |<--| ELEMENTS OF MATRIX
!| A12            |<--| ELEMENTS OF MATRIX
!| A13            |<--| ELEMENTS OF MATRIX
!| A21            |<--| ELEMENTS OF MATRIX
!| A22            |<--| ELEMENTS OF MATRIX
!| A23            |<--| ELEMENTS OF MATRIX
!| A31            |<--| ELEMENTS OF MATRIX
!| A32            |<--| ELEMENTS OF MATRIX
!| A33            |<--| ELEMENTS OF MATRIX
!| F              |-->| FUNCTION USED IN THE FORMULA
!| ICOORD         |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
!| IKLE1          |-->| FIRST POINTS OF TRIANGLES
!| IKLE2          |-->| SECOND POINTS OF TRIANGLES
!| IKLE3          |-->| THIRD POINTS OF TRIANGLES
!| IKLE4          |-->| QUASI-BUBBLE POINT
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| SF             |-->| STRUCTURE OF FUNCTIONS F
!| SU             |-->| BIEF_OBJ STRUCTURE OF U
!| SURDET         |-->| HERE = 1/(2*SURFAC)
!| SV             |-->| BIEF_OBJ STRUCTURE OF V
!| U              |-->| FUNCTION U USED IN THE FORMULA
!| V              |-->| FUNCTION V USED IN THE FORMULA
!| XEL            |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
!| YEL            |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
!| XMUL           |-->| MULTIPLICATION FACTOR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_mt12aa => mt12aa
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
      INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
!
      DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
      DOUBLE PRECISION, INTENT(INOUT) :: A21(*),A22(*),A23(*)
      DOUBLE PRECISION, INTENT(INOUT) :: A31(*),A32(*),A33(*)
!
      DOUBLE PRECISION, INTENT(IN) :: XMUL
      DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
!
!     STRUCTURES OF F, U, V
      TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
!
      DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
      DOUBLE PRECISION, INTENT(IN) :: SURDET(nelmax)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,IELMF,IELMU,IELMV
      DOUBLE PRECISION X2,X3,Y2,Y3,F1,F2,F3,DEN
      DOUBLE PRECISION U1,U2,U3,V1,V2,V3,U123,V123,XSUR06,XSUR24
!
!-----------------------------------------------------------------------
!
      xsur06 = xmul/06.d0
      xsur24 = xmul/24.d0
!
!-----------------------------------------------------------------------
!
      ielmf=sf%ELM
      ielmu=su%ELM
      ielmv=sv%ELM
!
!  CASE WHERE F IS OF TYPE P1 AND U P0
!
      IF(ielmf.EQ.11.AND.ielmu.EQ.10.AND.ielmv.EQ.10) THEN
!
!================================
!  DERIVATIVE WRT X  =
!================================
!
      IF(icoord.EQ.1) THEN
!
!   LOOP ON THE ELEMENTS
!
      DO ielem = 1 , nelem
!
!   INITIALISES THE GEOMETRICAL VARIABLES
!
      x2 = xel(ielem,2)
      x3 = xel(ielem,3)
      y2 = yel(ielem,2)
      y3 = yel(ielem,3)
!
      f1  =  f(ikle1(ielem))
      f2  =  f(ikle2(ielem)) - f1
      f3  =  f(ikle3(ielem)) - f1
!
      den = (f3*y2 - f2*y3) * xsur06 * surdet(ielem)
!
!   EXTRADIAGONAL TERMS
!
      a23(ielem) = (  x3 *v(ielem) -  y3    *u(ielem) )*den
      a31(ielem) =-(  x2 *v(ielem) -     y2 *u(ielem) )*den
!
      a12(ielem) = - a23(ielem) - a31(ielem)
      a13(ielem) =   a12(ielem)
      a21(ielem) =   a23(ielem)
      a32(ielem) =   a31(ielem)
!
! DIAGONAL TERMS (THE SUM OF EACH COLUMN IS 0)
!
      a11(ielem) = - a21(ielem) - a31(ielem)
      a22(ielem) = - a12(ielem) - a32(ielem)
      a33(ielem) = - a13(ielem) - a23(ielem)
!
      ENDDO ! IELEM
!
      ELSEIF(icoord.EQ.2) THEN
!
!================================
!  DERIVATIVE WRT Y  =
!================================
!
      DO ielem = 1 , nelem
!
!   INITIALISES THE GEOMETRICAL VARIABLES
!
      x2  =  xel(ielem,2)
      x3  =  xel(ielem,3)
      y2  =  yel(ielem,2)
      y3  =  yel(ielem,3)
!
      f1  =  f(ikle1(ielem))
      f2  =  f(ikle2(ielem)) - f1
      f3  =  f(ikle3(ielem)) - f1
!
      den = (f3*x2 - f2*x3) * xsur06 * surdet(ielem)
!
!   EXTRADIAGONAL TERMS
!
      a23(ielem) = -( x3*v(ielem) - y3*u(ielem) ) * den
      a31(ielem) =  ( x2*v(ielem) - y2*u(ielem) ) * den
!
      a12(ielem) = - a23(ielem) - a31(ielem)
      a13(ielem) =   a12(ielem)
      a21(ielem) =   a23(ielem)
      a32(ielem) =   a31(ielem)
!
! DIAGONAL TERMS (THE SUM OF EACH COLUMN IS 0)
!
      a11(ielem) = - a21(ielem) - a31(ielem)
      a22(ielem) = - a12(ielem) - a32(ielem)
      a33(ielem) = - a13(ielem) - a23(ielem)
!
      ENDDO ! IELEM
!
        ELSE
!
          WRITE(lu,201) icoord
201       FORMAT(1x,'MT12AA (BIEF) : IMPOSSIBLE COMPONENT ',
     &              1i6,' CHECK ICOORD')
          CALL plante(0)
          stop
!
        ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(ielmf.EQ.11.AND.ielmu.EQ.11) THEN
!
!================================
!  DERIVATIVE WRT X  =
!================================
!
      IF(icoord.EQ.1) THEN
!
!   LOOP ON THE ELEMENTS
!
      DO ielem = 1 , nelem
!
!   INITIALISES THE GEOMETRICAL VARIABLES
!
      x2 = xel(ielem,2)
      x3 = xel(ielem,3)
      y2 = yel(ielem,2)
      y3 = yel(ielem,3)
!
      f1  =  f(ikle1(ielem))
      f2  =  f(ikle2(ielem)) - f1
      f3  =  f(ikle3(ielem)) - f1
!
      u1  =  u(ikle1(ielem))
      u2  =  u(ikle2(ielem))
      u3  =  u(ikle3(ielem))
!
      v1  =  v(ikle1(ielem))
      v2  =  v(ikle2(ielem))
      v3  =  v(ikle3(ielem))
!
      u123 = u1 + u2 + u3
      v123 = v1 + v2 + v3
!
      den = (f3*y2 - f2*y3) * xsur24 * surdet(ielem)
!
!   EXTRADIAGONAL TERMS
!
      a12(ielem) = ( (x2-x3)*(v123+v2) + (y3-y2)*(u123+u2) )*den
!
      a13(ielem) = ( (x2-x3)*(v123+v3) + (y3-y2)*(u123+u3) )*den
!
      a23(ielem) = ( x3*(v123+v3) - y3*(u123+u3) )*den
!
      a21(ielem) = ( x3*(v123+v1) - y3*(u123+u1) )*den
!
      a31(ielem) =-( x2*(v123+v1) - y2*(u123+u1) )*den
!
      a32(ielem) =-( x2*(v123+v2) - y2*(u123+u2) )*den
!
! DIAGONAL TERMS (THE SUM OF EACH COLUMN IS 0)
!
      a11(ielem) = - a21(ielem) - a31(ielem)
      a22(ielem) = - a12(ielem) - a32(ielem)
      a33(ielem) = - a13(ielem) - a23(ielem)
!
      ENDDO ! IELEM
!
      ELSEIF(icoord.EQ.2) THEN
!
!================================
!  DERIVATIVE WRT Y  =
!================================
!
      DO ielem = 1 , nelem
!
!   INITIALISES THE GEOMETRICAL VARIABLES
!
      x2  =  xel(ielem,2)
      x3  =  xel(ielem,3)
      y2  =  yel(ielem,2)
      y3  =  yel(ielem,3)
!
      f1  =  f(ikle1(ielem))
      f2  =  f(ikle2(ielem)) - f1
      f3  =  f(ikle3(ielem)) - f1
!
      u1  =  u(ikle1(ielem))
      u2  =  u(ikle2(ielem))
      u3  =  u(ikle3(ielem))
!
      v1  =  v(ikle1(ielem))
      v2  =  v(ikle2(ielem))
      v3  =  v(ikle3(ielem))
!
      u123 = u1 + u2 + u3
      v123 = v1 + v2 + v3
!
      den = (f3*x2 - f2*x3) * xsur24 * surdet(ielem)
!
!   EXTRADIAGONAL TERMS
!
      a12(ielem) =-( (x2-x3)*(v123+v2) + (y3-y2)*(u123+u2) )*den
!
      a13(ielem) =-( (x2-x3)*(v123+v3) + (y3-y2)*(u123+u3) )*den
!
      a23(ielem) =-( x3*(v123+v3) - y3*(u123+u3) )*den
!
      a21(ielem) =-( x3*(v123+v1) - y3*(u123+u1) )*den
!
      a31(ielem) = ( x2*(v123+v1) - y2*(u123+u1) )*den
!
      a32(ielem) = ( x2*(v123+v2) - y2*(u123+u2) )*den
!
! DIAGONAL TERMS (THE SUM OF EACH COLUMN IS 0)
!
      a11(ielem) = - a21(ielem) - a31(ielem)
      a22(ielem) = - a12(ielem) - a32(ielem)
      a33(ielem) = - a13(ielem) - a23(ielem)
!
      ENDDO ! IELEM
!
        ELSE
          WRITE(lu,201) icoord
          CALL plante(0)
          stop
        ENDIF
!
!     OTHER TYPES OF F FUNCTIONS
!
!-----------------------------------------------------------------------
!
      ELSE
        WRITE(lu,101) ielmf,ielmu
101     FORMAT(1x,'MT12AA (BIEF) :',/,
     &        1x,'COMBINATION OF F AND U: ',1i6,2x,1i6,' NOT AVAILABLE')
        CALL plante(0)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END