vc13aa.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE vc13aa
!                   *****************
!
     &( xmul,sf,f,xel,yel,
     &  ikle1,ikle2,ikle3,ikle4,nelem,nelmax,w1,w2,w3 , icoord )
!
!***********************************************************************
! BIEF   V7P0                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
!code
!+   (EXAMPLE OF THE X COMPONENT, WHICH CORRESPONDS TO ICOORD=1)
!+
!+                       /            DF
!+    VEC(I)  =  XMUL   /     ( P  *( --  )) D(OMEGA)
!+                     /OMEGA    I    DX
!+
!+    P   IS A LINEAR BASE
!+     I
!+
!+    F IS A VECTOR OF TYPE P1 OR OTHER
!
!note     IMPORTANT : IF F IS OF TYPE P0, THE RESULT IS 0.
!+                     HERE, IF F IS P0, IT REALLY MEANS THAT F IS
!+                     P1, BUT GIVEN BY ELEMENTS.
!+                     THE SIZE OF F SHOULD THEN BE : F(NELMAX,3).
!
!warning  THE JACOBIAN MUST BE POSITIVE
!warning  THE RESULT IS IN W IN NOT ASSEMBLED FORM
!
!history  J-M HERVOUET (LNH)    ; F LEPEINTRE (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 (EDF LAB, LNHE)
!+        12/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.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| F              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| ICOORD         |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
!| IKLE1          |-->| FIRST POINT OF TRIANGLES
!| IKLE2          |-->| SECOND POINT OF TRIANGLES
!| IKLE3          |-->| THIRD POINT OF TRIANGLES
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| SF             |-->| BIEF_OBJ STRUCTURE OF F
!| SURFAC         |-->| AREA OF TRIANGLES
!| W1             |<--| RESULT IN NON ASSEMBLED FORM
!| W2             |<--| RESULT IN NON ASSEMBLED FORM
!| W3             |<--| RESULT IN NON ASSEMBLED FORM
!| XEL            |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
!| XMUL           |-->| MULTIPLICATION COEFFICIENT
!| YEL            |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_vc13aa => vc13aa
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
      INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
      INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
!
      DOUBLE PRECISION, INTENT(IN)    :: XEL(nelmax,*),YEL(nelmax,*)
      DOUBLE PRECISION, INTENT(INOUT) ::W1(nelmax),W2(nelmax),W3(nelmax)
      DOUBLE PRECISION, INTENT(IN)    :: XMUL
!
!     STRUCTURE OF F AND REAL DATA
!
      TYPE(bief_obj), INTENT(IN) :: SF
      DOUBLE PRECISION, INTENT(IN) :: F(*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,IELMF
      DOUBLE PRECISION XSUR6,XSUR18,F1,F2,F3,F4,X2,X3,Y2,Y3
!
!-----------------------------------------------------------------------
!
      xsur6 = xmul / 6.d0
      xsur18= xmul /18.d0
!
!-----------------------------------------------------------------------
!
      ielmf=sf%ELM
!
!-----------------------------------------------------------------------
!
!     F IS LINEAR
!
      IF(ielmf.EQ.11) THEN
!
!     X COORDINATE
!
      IF(icoord.EQ.1) THEN
!
      DO ielem = 1 , nelem
!
        w1(ielem) = ( yel(ielem,2) *
     &                (f(ikle1(ielem))-f(ikle3(ielem)))
     &              + yel(ielem,3) *
     &                (f(ikle2(ielem))-f(ikle1(ielem))) ) * xsur6
        w2(ielem) = w1(ielem)
        w3(ielem) = w1(ielem)
!
      ENDDO
!
      ELSEIF(icoord.EQ.2) THEN
!
!     Y COORDINATE
!
      DO ielem = 1 , nelem
!
        w1(ielem) = ( xel(ielem,2) *
     &                (f(ikle3(ielem))-f(ikle1(ielem)))
     &              + xel(ielem,3) *
     &                (f(ikle1(ielem))-f(ikle2(ielem))) ) * xsur6
        w2(ielem)  =  w1(ielem)
        w3(ielem)  =  w1(ielem)
!
      ENDDO
!
      ELSE
!
        WRITE(lu,201) icoord
201     FORMAT(1x,'VC13AA (BIEF) : IMPOSSIBLE COMPONENT ',
     &            1i6,' CHECK ICOORD')
        CALL plante(0)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!     F IS QUASI-BUBBLE
!
      ELSEIF(ielmf.EQ.12) THEN
!
!     X COORDINATE
!
      IF(icoord.EQ.1) THEN
!
      DO ielem = 1 , nelem
!
        y2 = yel(ielem,2)
        y3 = yel(ielem,3)
!
        f1 = f(ikle1(ielem))
        f2 = f(ikle2(ielem)) - f1
        f3 = f(ikle3(ielem)) - f1
        f4 = f(ikle4(ielem)) - f1
!
        w1(ielem)=(y2*(-3*f4-2*f3+f2)+y3*(3*f4-f3+2*f2)) * xsur18
        w2(ielem)=(-3*y2*f3+y3*(-3*f4+f3+4*f2)) * xsur18
        w3(ielem)=(y2*(3*f4-4*f3-f2)+3*y3*f2) * xsur18
!
      ENDDO
!
      ELSEIF(icoord.EQ.2) THEN
!
!     Y COORDINATE
!
      DO ielem = 1 , nelem
!
        x2 = xel(ielem,2)
        x3 = xel(ielem,3)
!
        f1 = f(ikle1(ielem))
        f2 = f(ikle2(ielem)) - f1
        f3 = f(ikle3(ielem)) - f1
        f4 = f(ikle4(ielem)) - f1
!
        w1(ielem)=(x2*(3*f4+2*f3-f2)+x3*(-3*f4+f3-2*f2)) * xsur18
        w2(ielem)=(3*x2*f3+x3*(3*f4-f3-4*f2)) * xsur18
        w3(ielem)=(x2*(-3*f4+4*f3+f2)-3*x3*f2) * xsur18
!
      ENDDO
!
      ELSE
!
        WRITE(lu,201) icoord
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!     BEWARE: HERE F IS QUASI-BUBBLE BUT DISCONTINUOUS BETWEEN THE ELEMENTS
!
      ELSEIF(ielmf.EQ.16) THEN
!
!     X COORDINATE
!
      IF(icoord.EQ.1) THEN
!
      DO ielem = 1 , nelem
!
        y2 = yel(ielem,2)
        y3 = yel(ielem,3)
!
        f1 = f(ielem)
        f2 = f(ielem+  nelmax)-f1
        f3 = f(ielem+2*nelmax)-f1
        f4 = f(ielem+3*nelmax)-f1
!
        w1(ielem)=(y2*(-3*f4-2*f3+f2)+y3*(3*f4-f3+2*f2)) * xsur18
        w2(ielem)=(-3*y2*f3+y3*(-3*f4+f3+4*f2)) * xsur18
        w3(ielem)=(y2*(3*f4-4*f3-f2)+3*y3*f2) * xsur18
!
      ENDDO
!
      ELSEIF(icoord.EQ.2) THEN
!
!     Y COORDINATE
!
      DO ielem = 1 , nelem
!
        x2 = xel(ielem,2)
        x3 = xel(ielem,3)
!
        f1 = f(ielem)
        f2 = f(ielem+  nelmax)-f1
        f3 = f(ielem+2*nelmax)-f1
        f4 = f(ielem+3*nelmax)-f1
!
        w1(ielem)=(x2*(3*f4+2*f3-f2)+x3*(-3*f4+f3-2*f2)) * xsur18
        w2(ielem)=(3*x2*f3+x3*(3*f4-f3-4*f2)) * xsur18
        w3(ielem)=(x2*(-3*f4+4*f3+f2)-3*x3*f2) * xsur18
!
      ENDDO
!
      ELSE
!
        WRITE(lu,201) icoord
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      ELSEIF(ielmf.EQ.15) THEN
!
!     BEWARE: HERE F IS LINEAR BUT DISCONTINUOUS BETWEEN THE ELEMENTS
!
!     X COORDINATE
!
      IF(icoord.EQ.1) THEN
!
      DO ielem = 1 , nelem
!
        w1(ielem) = ( yel(ielem,2) * (f(ielem)-f(ielem+2*nelmax))
     &              + yel(ielem,3) * (f(ielem+nelmax)-f(ielem)))* xsur6
        w2(ielem) = w1(ielem)
        w3(ielem) = w1(ielem)
!
      ENDDO
!
      ELSEIF(icoord.EQ.2) THEN
!
!     Y COORDINATE
!
      DO ielem = 1 , nelem
!
        w1(ielem) = ( xel(ielem,2) * (f(ielem+2*nelmax)-f(ielem))
     &              + xel(ielem,3) * (f(ielem)-f(ielem+nelmax)))*xsur6
        w2(ielem) =  w1(ielem)
        w3(ielem) =  w1(ielem)
!
      ENDDO
!
      ELSE
!
        WRITE(lu,201) icoord
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      ELSE
!
!-----------------------------------------------------------------------
!
        WRITE(lu,102) ielmf,sf%NAME
102     FORMAT(1x,'VC13AA (BIEF) :',/,
     &         1x,'DISCRETISATION OF F : ',1i6,' NOT IMPLEMENTED',/,
     &         1x,'REAL NAME OF F: ',a6)
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END