vc01tt.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE vc01tt
!                   *****************
!
     &( xmul,sf,f,x,y,z,
     &  ikle1,ikle2,ikle3,ikle4,nelem,nelmax,w1,w2,w3,w4)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
!code
!+                    /
!+    VEC(I) = XMUL  /    PSI(I) * F  D(OMEGA)
!+                  /OMEGA
!+
!+    PSI(I) IS A BASE OF TYPE P1 TETRAHEDRON
!+
!+    F IS A VECTOR OF TYPE IELMF
!
!warning  THE JACOBIAN MUST BE POSITIVE
!warning  THE RESULT IS IN W IN NOT ASSEMBLED FORM - REAL MESH
!
!history  J-M HERVOUET (LNH)
!+        22/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 VECTOR FORMULA
!| IKLE1          |-->| FIRST POINT OF TETRAHEDRA
!| IKLE2          |-->| SECOND POINT OF TETRAHEDRA
!| IKLE3          |-->| THIRD POINT OF TETRAHEDRA
!| IKLE4          |-->| FOURTH POINT OF TETRAHEDRA
!| 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
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH
!| XMUL           |-->| MULTIPLICATION COEFFICIENT
!| Y              |-->| ORDINATES OF POINTS IN THE MESH
!| Z              |-->| ELEVATIONS OF POINTS IN THE MESH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_vc01tt => vc01tt
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX
      INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax)
      INTEGER, INTENT(IN) :: IKLE3(nelmax),IKLE4(nelmax)
!
      DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),Z(*)
      DOUBLE PRECISION, INTENT(INOUT) :: W1(nelmax),W2(nelmax)
      DOUBLE PRECISION, INTENT(INOUT) :: W3(nelmax),W4(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 XSUR120,COEF,F1234,X2,X3,X4,Y2,Y3,Y4,Z2,Z3,Z4
      DOUBLE PRECISION F1,F2,F3,F4
      INTEGER I1,I2,I3,I4,IELEM,IELMF
!
!***********************************************************************
!
      ielmf=sf%ELM
!
!-----------------------------------------------------------------------
!
!   F IS LINEAR
!
      IF(ielmf.EQ.31.OR.ielmf.EQ.51) THEN
!
        xsur120 = xmul / 120.d0
!
        DO ielem = 1 , nelem
!
          i1 = ikle1(ielem)
          i2 = ikle2(ielem)
          i3 = ikle3(ielem)
          i4 = ikle4(ielem)
!
          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)
!
          z2 = z(i2)-z(i1)
          z3 = z(i3)-z(i1)
          z4 = z(i4)-z(i1)
!
          coef = xsur120*
     &          (x2*y3*z4-x2*y4*z3-y2*x3*z4+y2*x4*z3+z2*x3*y4-z2*x4*y3)
!
          f1  = f(ikle1(ielem))
          f2  = f(ikle2(ielem))
          f3  = f(ikle3(ielem))
          f4  = f(ikle4(ielem))
!
          f1234 = f1+f2+f3+f4
          w1(ielem) = coef * (f1234+f1)
          w2(ielem) = coef * (f1234+f2)
          w3(ielem) = coef * (f1234+f3)
          w4(ielem) = coef * (f1234+f4)
!
        ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
!   F IS CONSTANT BY ELEMENT
!
!   SAME METHOD BUT F HAS THE SAME VALUE
!   FOR THE 4 POINTS OF THE TETRAHEDRON
!
      ELSE IF(ielmf.EQ.30.OR.ielmf.EQ.50) THEN
!
        xsur120 = xmul / 120.d0
!
        DO ielem = 1 , nelem
!
          i1 = ikle1(ielem)
          i2 = ikle2(ielem)
          i3 = ikle3(ielem)
          i4 = ikle4(ielem)
!
          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)
!
          z2 = z(i2)-z(i1)
          z3 = z(i3)-z(i1)
          z4 = z(i4)-z(i1)
!
          coef = xsur120*
     &          (x2*y3*z4-x2*y4*z3-y2*x3*z4+y2*x4*z3+z2*x3*y4-z2*x4*y3)
!
! COULD BE SIMPLIFIED BUT IS NICE TO KEEP
! THE SAME FORM AS ABOVE
!
          f1  = f(ielem)
          f2  = f1
          f3  = f1
          f4  = f1
!
          f1234 = f1+f2+f3+f4
          w1(ielem) = coef * (f1234+f1)
          w2(ielem) = coef * (f1234+f2)
          w3(ielem) = coef * (f1234+f3)
          w4(ielem) = coef * (f1234+f4)
!
        ENDDO ! IELEM
!
!-----------------------------------------------------------------------
      ELSE
!
        WRITE(lu,102) ielmf,sf%NAME
102     FORMAT(1x,'VC01TT (BIEF):',/,
     &         1x,'DISCRETISATION OF F : ',1i6,' NOT IMPLEMENTED',/,
     &         1x,'REAL NAME OF F: ',a6)
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE vc01tt