kspg11.f

Go to the documentation of this file.

!                 *****************
                  SUBROUTINE kspg11
!                 *****************
!
     &(kx,ky,xel,yel,u,v,ikle,nelem,nelmax,xmul)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES A VECTOR, USED BY THE METHOD:
!+                STREAMLINE UPWIND PETROV GALERKIN (SUPG)
!+                WITH AN OFF-CENTERING OF 1.
!code
!+                    DX   U
!+             KX = -----------
!+                   2 NORM(U)
!+
!+                    DY   V
!+             KY = -----------
!+                   2 NORM(U)
!
!code
!+  MEANING OF IELM :
!+
!+  TYPE OF ELEMENT      NUMBER OF POINTS      CODED IN THIS SUBROUTINE
!+
!+  11 : P1 TRIANGLE            3                       YES
!+  12 : P2 TRIANGLE            6
!+  21 : Q1 QUADRILATERAL       4                       YES
!+  41 : TELEMAC-3D PRISMS      6
!
!history  J-M HERVOUET (LNH)
!+        08/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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IKLE           |-->| CONNECTIVITY TABLE.
!| KX             |-->| FIRST COMPONENT OF RESULTING VECTOR
!| KY             |-->| SECOND COMPONENT OF RESULTING VECTOR
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| U              |-->| FIRST COMPONENT OF VELOCITY
!| V              |-->| SECOND COMPONENT OF VELOCITY
!| XMUL           |-->| MULTIPLICATION COEFICIENT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NELEM,NELMAX
      INTEGER, INTENT(IN)             :: IKLE(nelmax,*)
      DOUBLE PRECISION, INTENT(INOUT) :: KX(nelem),KY(nelem)
      DOUBLE PRECISION, INTENT(IN)    :: U(*),V(*),XMUL
      DOUBLE PRECISION, INTENT(IN)    :: XEL(nelmax,*),YEL(nelmax,*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,I1,I2,I3
!
      DOUBLE PRECISION UMOY,VMOY,H,SUNORM,X2,X3,Y2,Y3
      DOUBLE PRECISION SURFAC,GP1X,GP1Y,GP2X,GP2Y,GP3X,GP3Y
      DOUBLE PRECISION A1,A2,A3,H1,H2,H3,C1,C2,C3,UNORM,VNORM
!
      INTRINSIC max,sqrt
!
!-----------------------------------------------------------------------
!
      DO ielem = 1 , nelem
!
        i1 = ikle(ielem,1)
        i2 = ikle(ielem,2)
        i3 = ikle(ielem,3)
!
        x2 = xel(ielem,2)
        x3 = xel(ielem,3)
        y2 = yel(ielem,2)
        y3 = yel(ielem,3)
!
        gp1x = y2-y3
        gp1y = x3-x2
        sunorm = 1.d0 / sqrt(gp1x**2+gp1y**2)
        gp1x = gp1x * sunorm
        gp1y = gp1y * sunorm
!
        gp2x = y3
        gp2y =   -x3
        sunorm = 1.d0 / sqrt(gp2x**2+gp2y**2)
        gp2x = gp2x * sunorm
        gp2y = gp2y * sunorm
!
        gp3x =   -y2
        gp3y = x2
        sunorm = 1.d0 / sqrt(gp3x**2+gp3y**2)
        gp3x = gp3x * sunorm
        gp3y = gp3y * sunorm
!
        c3 = sqrt(  x2**2     +  y2**2 )
        c1 = sqrt( (x3-x2)**2 + (y3-y2)**2 )
        c2 = sqrt(  x3**2     +  y3**2 )
!
        surfac = 0.5d0 * (x2*y3 - x3*y2)
!
        h1 = 2*surfac/c1
        h2 = 2*surfac/c2
        h3 = 2*surfac/c3
!
        h = max(h1,h2,h3)
!
        umoy = u(i1) + u(i2) + u(i3)
        vmoy = v(i1) + v(i2) + v(i3)
!
        sunorm = 1.d0 / max( sqrt(umoy**2+vmoy**2) , 1.d-10 )
!
        unorm = umoy * sunorm
        vnorm = vmoy * sunorm
!
        a1 = gp1x * unorm + gp1y * vnorm
        a2 = gp2x * unorm + gp2y * vnorm
        a3 = gp3x * unorm + gp3y * vnorm
!
        IF(a1*h.GT.h1) h = h1
        IF(a2*h.GT.h2) h = h2
        IF(a3*h.GT.h3) h = h3
!
        kx(ielem) = 0.33333333d0 * xmul * h * unorm
        ky(ielem) = 0.33333333d0 * xmul * h * vnorm
!
      ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
      RETURN
      END