mt02aa_2.f

Go to the documentation of this file.

!                   *******************
                    SUBROUTINE mt02aa_2
!                   *******************
!
     &( a11 , a12 , a13 ,
     &        a22 , a23 ,
     &              a33 ,
     &  xmul,su,sv,
     &  xel,yel,surfac,nelem,nelmax)
!
!***********************************************************************
! BIEF   V7P0                                   21/08/2010
!***********************************************************************
!
!brief    BUILDS THE DIFFUSION TERM FOR ESTEL2D.
!
!warning  THE JACOBIAN MUST BE POSITIVE
!
!history  J-M HERVOUET (LNH)    ; F  LEPEINTRE (LNH)
!+        28/11/94
!+        V5P8
!+
!
!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 (EDF LAB, LNHE)
!+        13/05/2014
!+        V7P0
!+   Discontinuous elements better treated: new types 15, 16 and 17 for
!+   discontinuous linear, quasi-bubble, and quadratic, rather than
!+   using component DIMDISC=11, 12 or 13.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A11            |<--| ELEMENTS OF MATRIX
!| A12            |<--| ELEMENTS OF MATRIX
!| A13            |<--| ELEMENTS OF MATRIX
!| A22            |<--| ELEMENTS OF MATRIX
!| A23            |<--| ELEMENTS OF MATRIX
!| A33            |<--| ELEMENTS OF MATRIX
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| SU             |-->| BIEF_OBJ STRUCTURE OF U
!| SURFAC         |-->| AREA OF TRIANGLES
!| SV             |-->| BIEF_OBJ STRUCTURE OF V
!| XEL            |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
!| YEL            |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
!| XMUL           |-->| MULTIPLICATION FACTOR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_mt02aa_2 => mt02aa_2
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX
      DOUBLE PRECISION, INTENT(INOUT) :: A11(*),A12(*),A13(*)
      DOUBLE PRECISION, INTENT(INOUT) ::        A22(*),A23(*)
      DOUBLE PRECISION, INTENT(INOUT) ::               A33(*)
      DOUBLE PRECISION, INTENT(IN) :: XMUL
!     STRUCTURE OF U AND V
      TYPE(bief_obj)  , INTENT(IN) :: SU,SV
      DOUBLE PRECISION, INTENT(IN) :: XEL(nelmax,3),YEL(nelmax,3)
      DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
!
      INTEGER IELEM,IELMNU,IELMNV
!
      DOUBLE PRECISION X2,X3,Y2,Y3,KSAT1,KSAT2,KSAT3,SOM,XSUR12
!
!=======================================================================
!
!     EXTRACTS THE TYPE OF ELEMENT FOR VISCOSITY
!
      ielmnu = su%ELM
      ielmnv = sv%ELM
!
      xsur12 = xmul / 12.d0
!
!-----------------------------------------------------------------------
! TESTS THE TYPES OF U AND V
! U (KR) : P0 AND DIM 3 (BECAUSE DISCONTINUOUS P1) - V (KS) : P0 AND DIM 3
!-----------------------------------------------------------------------
!
      IF(ielmnu.EQ.15.AND.ielmnv.EQ.15) THEN
!
      DO ielem = 1 , nelem
!
! THE 3 TERMS OF MATRIX V (KS IS SYMMETRICAL)
!
        ksat1=sv%R(ielem)
        ksat2=sv%R(ielem+nelem)
        ksat3=sv%R(ielem+2*nelem)
!
!   INITIALISES THE GEOMETRICAL VARIABLES
!
        x2  =  xel(ielem,2)
        x3  =  xel(ielem,3)
!
        y2  =  yel(ielem,2)
        y3  =  yel(ielem,3)
!
!   INITIALISES THE INTERMEDIATE VARIABLES
!
        som = ( su%R(ielem+2*nelem)
     &      +   su%R(ielem+nelem)
     &      +   su%R(ielem) ) * xsur12 / surfac(ielem)
!
!  DIAGONAL TERMS
!
        a11(ielem) = (ksat1*y2**2-2*ksat1*y2*y3+ksat1*y3**2+ksat2*x2**2-
     &  2*ksat2*x2*x3+ksat2*x3**2-2*ksat3*y2*x2+2*ksat3*y2*x3+
     &  2*ksat3*x2*y3-2*ksat3*y3*x3)*som
!
        a22(ielem) = (ksat1*y3**2+ksat2*x3**2-2*ksat3*y3*x3)*som
!
        a33(ielem) = (ksat1*y2**2+ksat2*x2**2-2*ksat3*y2*x2)*som
!
!  EXTRADIAGONAL TERMS
!
        a12(ielem) = -(-ksat1*y2*y3+ksat1*y3**2-ksat2*x2*x3+ksat2*x3**2+
     &          ksat3*x2*y3-2*ksat3*y3*x3+ksat3*y2*x3)*som
!
        a13(ielem) = -(ksat1*y2**2-ksat1*y2*y3+ksat2*x2**2-ksat2*x2*x3-
     &          2*ksat3*y2*x2+ksat3*y2*x3+ksat3*x2*y3)*som
!
        a23(ielem) = (-ksat1*y2*y3-ksat2*x2*x3+ksat3*y2*x3+ksat3*x2*y3)*
     &          som
!
!   END OF THE LOOP ON THE ELEMENTS
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      ELSE
!
        WRITE(lu,11)
11      FORMAT(1x,'MT02AA_2 (BIEF) : TYPES NOT AVAILABLE')
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END