vc08pp.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE vc08pp
!                   *****************
!
     &( xmul,sf,su,sv,sw,f,u,v,w,x,y,z,
     &  ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,nelem,nelmax,
     &  w1,w2,w3,w4,w5,w6 )
!
!***********************************************************************
! BIEF   V6P3                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE FOLLOWING VECTOR IN FINITE ELEMENTS:
!code
!+                    /                  DF      DF
!+      V  =  XMUL   /       PSII  * ( U --  + V -- )   D(OMEGA)
!+       I          /OMEGA               DX      DY
!+
!+    PSI(I) IS A BASE OF TYPE P1 PRISM
!
!warning  THE JACOBIAN MUST BE POSITIVE
!warning  THE RESULT IS IN W IN NOT ASSEMBLED FORM - REAL MESH
!
!history  J-M HERVOUET (LNH)
!+        23/02/96
!+        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
!| 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
!| SW             |-->| BIEF_OBJ STRUCTURE OF W
!| U              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| V              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| W              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| 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
!| Y              |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
!| XMUL           |-->| MULTIPLICATION COEFFICIENT
!| Z              |-->| ELEVATIONS OF POINTS ,PER POINT !!!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_vc08pp => vc08pp
!
      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),Z(*)
      DOUBLE PRECISION, INTENT(IN) :: XMUL
      DOUBLE PRECISION, INTENT(INOUT)::W1(nelmax),W2(nelmax),W3(nelmax)
      DOUBLE PRECISION, INTENT(INOUT)::W4(nelmax),W5(nelmax),W6(nelmax)
!
!     STRUCTURES OF F, U, V AND REAL DATA
!
      TYPE(bief_obj), INTENT(IN) :: SF,SU,SV,SW
      DOUBLE PRECISION, INTENT(IN) ::F(*),U(*),V(*),W(*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION X2,X3,Y2,Y3,DEN
      DOUBLE PRECISION PZ1,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 Q1,Q2,Q3,Q4,Q5,Q6,H1,H2,H3,HT
      DOUBLE PRECISION F1,F2,F3,F4,F5,F6
      DOUBLE PRECISION INT1,INT2,INT3,INT4,INT5,INT6
      DOUBLE PRECISION SUS,SUI,SHUI,SHUS
      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,SHVI,SHVS
      DOUBLE PRECISION SHVI1,SHVS1,SHVI2,SHVS2,SHVI3,SHVS3
      DOUBLE PRECISION HV1S3I,HV1I3S,HV1SI,HV2S3I,HV2I3S,HV2SI
      DOUBLE PRECISION HV3S3I,HV3I3S,HV3SI
      DOUBLE PRECISION W14,W41,W25,W52,W63,W36
!
      INTEGER I1,I2,I3,I4,I5,I6,IELEM,IELMF,IELMU,IELMV,IELMW
!
!**********************************************************************
!
      den = xmul/1440.d0
!
      ielmf=sf%ELM
      ielmu=su%ELM
      ielmv=sv%ELM
      ielmw=sw%ELM
!
!-----------------------------------------------------------------------
!
!     FUNCTION F AND VECTOR U ARE LINEAR
!
      IF(ielmu.EQ.41.AND.ielmv.EQ.41.AND.ielmw.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)
          q1  =  w(i1)
          q2  =  w(i2)
          q3  =  w(i3)
          q4  =  w(i4)
          q5  =  w(i5)
          q6  =  w(i6)
          f1  =  f(i1)
          f2  =  f(i2)
          f3  =  f(i3)
          f4  =  f(i4)
          f5  =  f(i5)
          f6  =  f(i6)
!
          h1  =  z(i4)-z(i1)
          h2  =  z(i5)-z(i2)
          h3  =  z(i6)-z(i3)
!
!    INTERMEDIATE COMPUTATIONS
!
          px2=y3*den
          px3=-y2*den
          px1=-px2-px3
          py2=-x3*den
          py3=x2*den
          py1=-py3-py2
          pz1=-den*(x2*y3-y2*x3)*2
!
          ht = h1+h2+h3
          sus = u6+u5+u4
          sui = u3+u2+u1
          shui = h3*u3+h2*u2+h1*u1
          shus = h3*u6+h2*u5+h1*u4
          shus1 = (ht+h1)*(sus+u4)+shus+u4*h1
          shui1 = (ht+h1)*(sui+u1)+shui+u1*h1
          shus2 = (ht+h2)*(sus+u5)+shus+u5*h2
          shui2 = (ht+h2)*(sui+u2)+shui+u2*h2
          shus3 = (ht+h3)*(sus+u6)+shus+u6*h3
          shui3 = (ht+h3)*(sui+u3)+shui+u3*h3
          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
          shvi = h3*v3+h2*v2+h1*v1
          shvs = h3*v6+h2*v5+h1*v4
          shvs1 = (ht+h1)*(svs+v4)+shvs+v4*h1
          shvi1 = (ht+h1)*(svi+v1)+shvi+v1*h1
          shvs2 = (ht+h2)*(svs+v5)+shvs+v5*h2
          shvi2 = (ht+h2)*(svi+v2)+shvi+v2*h2
          shvs3 = (ht+h3)*(svs+v6)+shvs+v6*h3
          shvi3 = (ht+h3)*(svi+v3)+shvi+v3*h3
          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
!
          w14 = q1+2*q4
          w41 = q4+2*q1
          w25 = q2+2*q5
          w52 = q5+2*q2
          w63 = q6+2*q3
          w36 = q3+2*q6
!
          int1=2*(3*w41+w52+w63)
          int2=2*(w41+w52)+w63
          int3=2*(w41+w63)+w52
          int4=2*(w41+3*w52+w63)
          int5=2*(w52+w63)+w41
          int6=2*(w41+w52+3*w63)
          w1(ielem)=pz1*((f1-f4)*int1+(f2-f5)*int2+(f3-f6)*int3)
          w2(ielem)=pz1*((f1-f4)*int2+(f2-f5)*int4+(f3-f6)*int5)
          w3(ielem)=pz1*((f1-f4)*int3+(f2-f4)*int5+(f3-f6)*int6)
!
          int1=2*(3*w14+w25+w36)
          int2=2*(w14+w25)+w36
          int3=2*(w14+w36)+w25
          int4=2*(w14+3*w25+w36)
          int5=2*(w25+w36)+w14
          int6=2*(w14+w25+3*w36)
          w4(ielem)=pz1*((f1-f4)*int1+(f2-f5)*int2+(f3-f6)*int3)
          w5(ielem)=pz1*((f1-f4)*int2+(f2-f5)*int4+(f3-f6)*int5)
          w6(ielem)=pz1*((f1-f4)*int3+(f2-f4)*int5+(f3-f6)*int6)
!
          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)=w1(ielem)+pxfi*hu1s3i+pyfi*hv1s3i
     &                       +pxfs* hu1si+pyfs* hv1si
          w2(ielem)=w2(ielem)+pxfi*hu2s3i+pyfi*hv2s3i
     &                       +pxfs* hu2si+pyfs* hv2si
          w3(ielem)=w3(ielem)+pxfi*hu3s3i+pyfi*hv3s3i
     &                       +pxfs* hu3si+pyfs* hv3si
          w4(ielem)=w4(ielem)+pxfi* hu1si+pyfi* hv1si
     &                       +pxfs*hu1i3s+pyfs*hv1i3s
          w5(ielem)=w5(ielem)+pxfi* hu2si+pyfi* hv2si
     &                       +pxfs*hu2i3s+pyfs*hv2i3s
          w6(ielem)=w6(ielem)+pxfi* hu3si+pyfi* hv3si
     &                       +pxfs*hu3i3s+pyfs*hv3i3s
!
        ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
      ELSE
!
!-----------------------------------------------------------------------
!
        WRITE(lu,101) ielmf,sf%NAME
        WRITE(lu,201) ielmu,su%NAME
        WRITE(lu,301)
101     FORMAT(1x,'VC08PP (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