vc18pp.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE vc18pp
!                   *****************
!
     &(xmul,sf,su,sv,f,u,v,x,y,
     & ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,nelem,nelmax,w1,w2,w3)
!
!***********************************************************************
! BIEF   V6P3                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE FOLLOWING VECTOR:
!code
!+                        /   1                  DF      DF
!+      V      =  XMUL   /   ---  PSI(IH)  * ( U --  + V -- )   D(OMEGA)
!+       IH,P*          /P*  DZ*                 DX      DY
!+
!+
!+    PSI(I) IS A 2D P1 BASE ON THE BASE TRIANGLE OF THE PRISM!
!+
!+    F AND U ARE 3D LINEAR BY PRISM
!
!warning  THE JACOBIAN MUST BE POSITIVE
!warning  THE RESULT IS IN W IN NOT ASSEMBLED FORM
!warning  COMPUTATION IN THE TRANSFORMED MESH !
!
!history  A. DECOENE (INRIA-LNHE)
!+        15/11/04
!+        V5P5
!+
!
!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
!| 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
!| SU             |-->| BIEF_OBJ STRUCTURE OF U
!| SV             |-->| BIEF_OBJ STRUCTURE OF V
!| SURFAC         |-->| AREA OF TRIANGLES
!| U              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| V              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| W1             |<--| RESULT IN NON ASSEMBLED FORM
!| W2             |<--| RESULT IN NON ASSEMBLED FORM
!| W3             |<--| RESULT IN NON ASSEMBLED FORM
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
!| Y              |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
!| XMUL           |-->| MULTIPLICATION COEFFICIENT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_vc18pp => vc18pp
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX
      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),XMUL
      DOUBLE PRECISION, INTENT(INOUT) ::W1(nelmax),W2(nelmax),W3(nelmax)
!
!     STRUCTURES OF F, U, V AND REAL DATA
!
      TYPE(bief_obj), INTENT(IN) :: SF,SU,SV
      DOUBLE PRECISION, INTENT(IN) :: F(*),U(*),V(*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,IELMF,IELMU,IELMV
      DOUBLE PRECISION X2,X3,Y2,Y3,DEN
      DOUBLE PRECISION PX1,PX2,PX3,PY1,PY2,PY3
      DOUBLE PRECISION PXFI,PXFS,PYFI,PYFS
      DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6
      DOUBLE PRECISION F1,F2,F3,F4,F5,F6
      DOUBLE PRECISION SUS,SUI
      DOUBLE PRECISION SHUI1,SHUS1,SHUI2,SHUS2,SHUI3,SHUS3
      DOUBLE PRECISION HU1S3I,HU1I3S,HU1SI,HU2S3I,HU2I3S,HU2SI
      DOUBLE PRECISION HU3S3I,HU3I3S,HU3SI
      DOUBLE PRECISION SVS,SVI
      DOUBLE PRECISION SHVI1,SHVS1,SHVI2,SHVS2,SHVI3,SHVS3
      DOUBLE PRECISION HV1S3I,HV1I3S,HV1SI,HV2S3I,HV2I3S,HV2SI
      DOUBLE PRECISION HV3S3I,HV3I3S,HV3SI
!
      INTEGER I1,I2,I3,I4,I5,I6
!
!**********************************************************************
!
      den = xmul/1440.d0
!
      ielmf=sf%ELM
      ielmu=su%ELM
      ielmv=sv%ELM
!
!-----------------------------------------------------------------------
!
!     FUNCTION F AND VECTOR U ARE LINEAR
!
      IF(ielmu.EQ.41.AND.ielmv.EQ.41.AND.ielmf.EQ.41) THEN
!
!        LOOP ON THE ELEMENTS
!
        DO ielem = 1,nelem
!
          i1 = ikle1(ielem)
          i2 = ikle2(ielem)
          i3 = ikle3(ielem)
          i4 = ikle4(ielem)
          i5 = ikle5(ielem)
          i6 = ikle6(ielem)
!
!         X2  =   X(I2) - X(I1)
!         X3  =   X(I3) - X(I1)
!         Y2  =   Y(I2) - Y(I1)
!         Y3  =   Y(I3) - Y(I1)
!
          x2 = x(ielem,2)
          x3 = x(ielem,3)
          y2 = y(ielem,2)
          y3 = y(ielem,3)
!
          u1  =  u(i1)
          u2  =  u(i2)
          u3  =  u(i3)
          u4  =  u(i4)
          u5  =  u(i5)
          u6  =  u(i6)
          v1  =  v(i1)
          v2  =  v(i2)
          v3  =  v(i3)
          v4  =  v(i4)
          v5  =  v(i5)
          v6  =  v(i6)
          f1  =  f(i1)
          f2  =  f(i2)
          f3  =  f(i3)
          f4  =  f(i4)
          f5  =  f(i5)
          f6  =  f(i6)
!
!    INTERMEDIATE COMPUTATIONS
!
          px2=y3*den
          px3=-y2*den
          px1=-px2-px3
          py2=-x3*den
          py3=x2*den
          py1=-py3-py2
!
          sus    = u6+u5+u4
          sui    = u3+u2+u1
          shus1  = 5 * (sus+u4)
          shui1  = 5 * (sui+u1)
          shus2  = 5 * (sus+u5)
          shui2  = 5 * (sui+u2)
          shus3  = 5 * (sus+u6)
          shui3  = 5 * (sui+u3)
          hu1s3i = shus1 + 3*shui1
          hu1i3s = shui1 + 3*shus1
          hu1si  = shus1 +   shui1
          hu2s3i = shus2 + 3*shui2
          hu2i3s = shui2 + 3*shus2
          hu2si  = shus2 +   shui2
          hu3s3i = shus3 + 3*shui3
          hu3i3s = shui3 + 3*shus3
          hu3si  = shus3 +   shui3
!
          svs    = v6+v5+v4
          svi    = v3+v2+v1
          shvs1  = 5 * (svs+v4)
          shvi1  = 5 * (svi+v1)
          shvs2  = 5 * (svs+v5)
          shvi2  = 5 * (svi+v2)
          shvs3  = 5 * (svs+v6)
          shvi3  = 5 * (svi+v3)
          hv1s3i = shvs1 + 3*shvi1
          hv1i3s = shvi1 + 3*shvs1
          hv1si  = shvs1 +   shvi1
          hv2s3i = shvs2 + 3*shvi2
          hv2i3s = shvi2 + 3*shvs2
          hv2si  = shvs2 +   shvi2
          hv3s3i = shvs3 + 3*shvi3
          hv3i3s = shvi3 + 3*shvs3
          hv3si  = shvs3 +   shvi3
!
          pxfi=px1*f1+px2*f2+px3*f3
          pxfs=px1*f4+px2*f5+px3*f6
          pyfi=py1*f1+py2*f2+py3*f3
          pyfs=py1*f4+py2*f5+py3*f6
!
          w1(ielem)= pxfi * (hu1s3i+hu1si)
     &             + pyfi * (hv1s3i+hv1si)
     &             + pxfs * (hu1si+hu1i3s)
     &             + pyfs * (hv1si+hv1i3s)
!
          w2(ielem)= pxfi * (hu2s3i+hu2si)
     &             + pyfi * (hv2s3i+hv2si)
     &             + pxfs * (hu2si+hu2i3s)
     &             + pyfs * (hv2si+hv2i3s)
!
          w3(ielem)= pxfi * (hu3s3i+hu3si)
     &             + pyfi * (hv3s3i+hv3si)
     &             + pxfs * (hu3si+hu3i3s)
     &             + pyfs * (hv3si+hv3i3s)
!
        ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
      ELSE
!
!-----------------------------------------------------------------------
!
        WRITE(lu,101) ielmf,sf%NAME
        WRITE(lu,201) ielmu,su%NAME
        WRITE(lu,301)
101     FORMAT(1x,'VC18PP (BIEF) :',/,
     &         1x,'DISCRETIZATION OF F:',1i6,
     &         1x,'REAL NAME: ',a6)
201     FORMAT(1x,'DISCRETIZATION OF U:',1i6,
     &         1x,'REAL NAME: ',a6)
301     FORMAT(1x,'CASE NOT IMPLEMENTED')
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END