vc13pp2.f

Go to the documentation of this file.

!                   ******************
                    SUBROUTINE vc13pp2
!                   ******************
!
     &(xmul,sf,f,x,y,z,
     & ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,nelem,nelmax,
     & w1,w2,w3,w4,w5,w6,icoord)
!
!***********************************************************************
! BIEF   V6P3                                   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
!+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!+    HERE F IS A FUNCTION OF X AND Y (NOT OF Z)
!+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!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 R&D LNHE)
!+        07/01/2013
!+        V6P3
!+   X and Y are now given per element.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| F              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| FORMUL         |-->| SEE AT THE END OF THE SUBROUTINE
!| ICOORD         |-->| 1: DERIVATIVE ALONG X, 2: ALONG Y
!| IKLE1          |-->| FIRST POINT OF PRISMS
!| IKLE2          |-->| SECOND POINT OF PRISMS
!| IKLE3          |-->| THIRD POINT OF PRISMS
!| IKLE4          |-->| FOURTH POINT OF PRISMS
!| IKLE5          |-->| FIFTH POINT OF PRISMS
!| IKLE6          |-->| SIXTH POINT OF PRISMS
!| 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
!| W4             |<--| RESULT IN NON ASSEMBLED FORM
!| W5             |<--| RESULT IN NON ASSEMBLED FORM
!| W6             |<--| RESULT IN NON ASSEMBLED FORM
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
!| XMUL           |-->| MULTIPLICATION COEFFICIENT
!| Y              |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
!| Z              |-->| Z-COORDINATES OF POINTS IN THE MESH, PER POINT !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_vc13pp2 => vc13pp2
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,ICOORD
      INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
      INTEGER, INTENT(IN) :: IKLE4(nelmax),IKLE5(nelmax),IKLE6(nelmax)
!
      DOUBLE PRECISION, INTENT(IN)    :: X(nelmax,6),Y(nelmax,6),Z(*)
      DOUBLE PRECISION, INTENT(INOUT) ::W1(nelmax),W2(nelmax),W3(nelmax)
      DOUBLE PRECISION, INTENT(INOUT) ::W4(nelmax),W5(nelmax),W6(nelmax)
      DOUBLE PRECISION, INTENT(IN)    :: XMUL
!
!     STRUCTURE OF F AND REAL DATA
!
      TYPE(bief_obj), INTENT(IN) :: SF
      DOUBLE PRECISION, INTENT(IN) :: F(*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION XSUR48,F1,F2,F3
      DOUBLE PRECISION X2,X3,Y2,Y3,H1,H2,H3,H123,FX6,FY6
      INTEGER I1,I2,I3,I4,I5,I6,IELEM,IELMF
!
!-----------------------------------------------------------------------
!
      xsur48  = xmul/48.d0
!
!-----------------------------------------------------------------------
!
      ielmf=sf%ELM
!
!=======================================================================
!
!     F IS LINEAR
!
      IF(ielmf.EQ.41) THEN
!
      IF(icoord.EQ.1) THEN
!
!-----------------------------------------------------------------------
!
!  DERIVATIVE WRT X
!
      DO ielem = 1 , nelem
!
        i1 = ikle1(ielem)
        i2 = ikle2(ielem)
        i3 = ikle3(ielem)
        i4 = ikle4(ielem)
        i5 = ikle5(ielem)
        i6 = ikle6(ielem)
!
        f1 = f(i1)
        f2 = f(i2)
        f3 = f(i3)
!
!       REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
!
        y2  =  y(ielem,2)
        y3  =  y(ielem,3)
!
        fx6 = y2*(f1-f3) + y3*(f2-f1)
!
        h1 = z(i4) - z(i1)
        h2 = z(i5) - z(i2)
        h3 = z(i6) - z(i3)
        h123 = h1 + h2 + h3
!
        w1(ielem) = xsur48 * (h123+h1) * fx6
        w2(ielem) = xsur48 * (h123+h2) * fx6
        w3(ielem) = xsur48 * (h123+h3) * fx6
!
        w4(ielem) = w1(ielem)
        w5(ielem) = w2(ielem)
        w6(ielem) = w3(ielem)
!
      ENDDO ! IELEM
!
      ELSEIF(icoord.EQ.2) THEN
!
!-----------------------------------------------------------------------
!
!  DERIVATIVE WRT Y
!
      DO ielem = 1 , nelem
!
        i1 = ikle1(ielem)
        i2 = ikle2(ielem)
        i3 = ikle3(ielem)
        i4 = ikle4(ielem)
        i5 = ikle5(ielem)
        i6 = ikle6(ielem)
!
        f1 = f(i1)
        f2 = f(i2)
        f3 = f(i3)
!
!       REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
!
        x2  =  x(ielem,2)
        x3  =  x(ielem,3)
!
        fy6 = x2*(f3-f1) + x3*(f1-f2)
!
        h1 = z(i4) - z(i1)
        h2 = z(i5) - z(i2)
        h3 = z(i6) - z(i3)
        h123 = h1 + h2 + h3
!
        w1(ielem) = xsur48 * (h123+h1) * fy6
        w2(ielem) = xsur48 * (h123+h2) * fy6
        w3(ielem) = xsur48 * (h123+h3) * fy6
!
        w4(ielem) = w1(ielem)
        w5(ielem) = w2(ielem)
        w6(ielem) = w3(ielem)
!
        ENDDO
!
      ELSE
!
!-----------------------------------------------------------------------
!
        WRITE(lu,201) icoord
201     FORMAT(1x,'VC13PP2 (BIEF) : IMPOSSIBLE COMPONENT ',
     &            1i6,' CHECK ICOORD')
        CALL plante(1)
        stop
!
      ENDIF
!
!=======================================================================
!
      ELSE
!
!=======================================================================
!
        WRITE(lu,102) ielmf,sf%NAME
102     FORMAT(1x,'VC13PP2 (BIEF) :',/,
     &         1x,'DISCRETISATION OF F : ',1i6,' NOT IMPLEMENTED',/,
     &         1x,'REAL NAME OF F: ',a6)
        CALL plante(1)
        stop
!
      ENDIF
!
!=======================================================================
!
      RETURN
      END