vc00tt.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE vc00tt
!                   *****************
!
     &(xmul,x,y,z,ikle1,ikle2,ikle3,ikle4,
     & nelem,nelmax,w1,w2,w3,w4,formul,npoin2,nelem2,ielm1)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
!code
!+                    /
!+    VEC(I) = XMUL  /    PSI(I)  D(OMEGA)
!+                  /OMEGA
!+
!+    PSI(I) IS A BASE OF TYPE P1 TETRAHEDRON
!
!warning  THE JACOBIAN MUST BE POSITIVE
!warning  THE RESULT IS IN W IN NOT ASSEMBLED FORM
!
!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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FORMUL         |-->| STRING WITH THE FORMULA DESCRIBING THE VECTOR
!| IELM1          |-->| TYPE OF ELEMENT
!| 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
!| NPOIN2         |-->| NUMBER OF POINTS IN THE UNDERLYING 2D MESH
!| 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 declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN2,NELEM2,IELM1
      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)
      DOUBLE PRECISION, INTENT(INOUT) :: W2(nelmax)
      DOUBLE PRECISION, INTENT(INOUT) :: W3(nelmax)
      DOUBLE PRECISION, INTENT(INOUT) :: W4(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XMUL
!
      CHARACTER(LEN=16), INTENT(IN) :: FORMUL
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION XSUR24,X2,X3,X4,Y2,Y3,Y4,Z2,Z3,Z4,H1,H2,H3,H4
      DOUBLE PRECISION COEF
      INTEGER I1,I2,I3,I4,IELEM,IP,I12D,I22D,I32D,I42D,IELEM2D
      INTEGER IT1,IT2,IT3
!
!-----------------------------------------------------------------------
!
      xsur24 = xmul/24.d0
!
!-----------------------------------------------------------------------
!
      IF(formul(1:7).EQ.'MASBAS ') THEN
!
!     STANDARD FORMULA
!
      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)
!
        w1(ielem) =
     &       (x2*(y3*z4-y4*z3)+y2*(x4*z3-x3*z4)+z2*(x3*y4-x4*y3))*xsur24
        w2(ielem) = w1(ielem)
        w3(ielem) = w1(ielem)
        w4(ielem) = w1(ielem)
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      ELSEIF(formul(1:7).EQ.'MASBAS2'.AND.ielm1.EQ.51) THEN
!
!     FORMULA WITH MASS-LUMPING
!
!     LOOP ON THE ELEMENTS
!
      DO ielem = 1 , nelem
!
        i1 = ikle1(ielem)
        i2 = ikle2(ielem)
        i3 = ikle3(ielem)
        i4 = ikle4(ielem)
!
!       RETRIEVING THE LOWER PLANE NUMBER
!
        ip=(min(i1,i2,i3,i4)-1)/npoin2 +1
!
!       RETRIEVING THE 2D POINTS NUMBERS ON THE SAME VERTICAL
!
        i12d=mod(i1-1,npoin2)+1
        i22d=mod(i2-1,npoin2)+1
        i32d=mod(i3-1,npoin2)+1
        i42d=mod(i4-1,npoin2)+1
!
!       RETRIEVING THE ORIGINAL PRISM HEIGHTS ON THE VERTICAL
!       TWO OF THE FOUR HEIGHTS WILL BE EQUAL, THE TWO CORRESPONDING
!       POINTS OF THE TETRAHEDRON BEING ON THE SAME VERTICAL
!
        h1=z(ip*npoin2+i12d)-z((ip-1)*npoin2+i12d)
        h2=z(ip*npoin2+i22d)-z((ip-1)*npoin2+i22d)
        h3=z(ip*npoin2+i32d)-z((ip-1)*npoin2+i32d)
        h4=z(ip*npoin2+i42d)-z((ip-1)*npoin2+i42d)
!
        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)
!
!       RETRIEVING THE CORRESPONDING TRIANGLE
        ielem2d=mod(ielem-1,nelem2)+1
!       IT1,IT2 AND IT3 FORM THE PROJECTION OF THE TETRAHEDRON ON THE
!       PLANE BELOW IT, IT IS A TRIANGLE
        it1=ikle1(ielem2d)+(ip-1)*npoin2
        it2=ikle2(ielem2d)+(ip-1)*npoin2
        it3=ikle3(ielem2d)+(ip-1)*npoin2
        x2 = x(it2)-x(it1)
        x3 = x(it3)-x(it1)
        y2 = y(it2)-y(it1)
        y3 = y(it3)-y(it1)
!
!       HERE COEF SHOULD BE PROPORTIONNAL TO THE VOLUME OF THE ELEMENT
!       BUT THE 3 TETRAHEDRONS IN A PRISM DO NOT HAVE THE SAME VOLUME
!       THEN THE SUM OF 3D VOLUMES ON THE VERTICAL DO NOT SUM CORRECTLY
!       TO GIVE DEPTH * VOLU2D. THE FORMULA BELOW ENSURES THIS PROPERTY
!
        coef=(x2*y3-x3*y2)*xsur24
!
        w1(ielem) = coef*h1
        w2(ielem) = coef*h2
        w3(ielem) = coef*h3
        w4(ielem) = coef*h4
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      ELSE
!
        WRITE(lu,*) ' '
        WRITE(lu,*) 'UNKNOWN FORMULA IN VC00TT:',formul
        WRITE(lu,*) 'WITH ELEMENT:',ielm1
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END