xref: /phasta/phSolver/compressible/e3q.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

        subroutine e3q (ycl,      shp,     shgl,
     &                  xl,       ql,      rmassl,
     &                  xmudmi,   sgn)
c
c----------------------------------------------------------------------
c
c This routine computes the element contribution to the
c diffusive flux vector and the lumped mass matrix.
c
c input:
c  ycl     (npro,nshl,ndof)       : Y variables
c  shp    (nen,ngauss)            : element shape-functions
c  shgl   (nsd,nen,ngauss)        : element local-grad-shape-functions
c  xl     (npro,nshl,nsd)        : nodal coordinates at current step
c
c output:
c  ql     (npro,nshl,idflx) : element RHS diffusion residual
c  rmassl     (npro,nshl)        : element lumped mass matrix
c
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension ycl(npro,nshl,ndof),
     &            shp(nshl,ngauss),
     &            shgl(nsd,nshl,ngauss),
     &            xl(npro,nenl,nsd),
     &            ql(npro,nshl,idflx),  rmassl(npro,nshl),
     &            xmudmi(npro,ngauss)
c
c local arrays
c
        dimension g1yi(npro,nflow),          g2yi(npro,nflow),
     &            g3yi(npro,nflow),          shg(npro,nshl,nsd),
     &            dxidx(npro,nsd,nsd),       WdetJ(npro),
     &            T(npro),                   cp(npro),
     &            u1(npro),                  u2(npro),
     &            u3(npro),                  rmu(npro),
     &            rlm(npro),                 rlm2mu(npro),
     &            con(npro),                 rho(npro)
c
        dimension qdi(npro,idflx),alph1(npro),alph2(npro)
c
        dimension sgn(npro,nshl),          shape(npro,nshl),
     &            shdrv(npro,nsd,nshl),    shpsum(npro)
c
c.... for surface tension
c
        dimension g1yti(npro),          g2yti(npro),
     &            g3yti(npro)
        integer idflow
c
c.... loop through the integration points
c

        ttim(33) = ttim(33) - secs(0.0)

        alph1 = 0.d0
        alph2 = 0.d0

        do intp = 1, ngauss
        if (Qwt(lcsyst,intp) .eq. zero) cycle          ! precaution
c
c.... create a matrix of shape functions (and derivatives) for each
c     element at this quadrature point. These arrays will contain
c     the correct signs for the hierarchic basis
c
        call getshp(shp,          shgl,      sgn,
     &              shape,        shdrv)
c
c.... initialize
c
        qdi = zero
c
c
c.... calculate the integration variables necessary for the
c     formation of q
c

        call e3qvar   (ycl,        shape,        shdrv,
     &                 rho,       xl,           g1yi,
     &                 g2yi,      g3yi,         shg,
     &                 dxidx,     WdetJ,        T,
     &                 cp,        u1,           u2,
     &                 u3                                 )
c
c.... compute diffusive flux vector at this integration point
c
c
c.... get material properties
c
        call getDiff (T,        cp,       rho,        ycl,
     &                rmu,      rlm,      rlm2mu,     con,  shape,
     &                xmudmi,   xl)

        idflow = 0
        if(idiff >= 1) then   !so taking care of all the idiff=1,2,3
        idflow = idflow+12
c
c.... compute diffusive fluxes
c
c.... diffusive flux in x1-direction
c
        qdi(:,1) =  rlm2mu      * g1yi(:,2)
     &               +      rlm * g2yi(:,3)
     &               +      rlm * g3yi(:,4)
        qdi(:,2) =  rmu         * g1yi(:,3)
     &               +      rmu * g2yi(:,2)
        qdi(:,3) =  rmu         * g1yi(:,4)
     &               +      rmu * g3yi(:,2)
        qdi(:,4) =  rlm2mu * u1 * g1yi(:,2) + rmu * u2 * g1yi(:,3)
     &                                   +    rmu * u3 * g1yi(:,4)
     &               +    rmu * u2 * g2yi(:,2) + rlm * u1 * g2yi(:,3)
     &               +    rmu * u3 * g3yi(:,2) + rlm * u1 * g3yi(:,4)
     &               +    con      * g1yi(:,5)

c
c.... diffusive flux in x2-direction
c
        qdi(:,5) =       rmu * g1yi(:,3)
     &            +      rmu * g2yi(:,2)
        qdi(:,6) =       rlm * g1yi(:,2)
     &            +   rlm2mu * g2yi(:,3)
     &            +      rlm * g3yi(:,4)
        qdi(:,7) =       rmu * g2yi(:,4)
     &            +      rmu * g3yi(:,3)
        qdi(:,8) =  rlm * u2 * g1yi(:,2) +    rmu * u1 * g1yi(:,3)
     &            + rmu * u1 * g2yi(:,2) + rlm2mu * u2 * g2yi(:,3)
     &            + rmu * u3 * g2yi(:,4)
     &            + rmu * u3 * g3yi(:,3) +    rlm * u2 * g3yi(:,4)
     &            +    con      * g2yi(:,5)
c
c.... diffusive flux in x3-direction
c
        qdi(:,9 ) =       rmu * g1yi(:,4)
     &            +      rmu * g3yi(:,2)
        qdi(:,10) =       rmu * g2yi(:,4)
     &            +      rmu * g3yi(:,3)
        qdi(:,11) =       rlm * g1yi(:,2)
     &            +      rlm * g2yi(:,3)
     &            +   rlm2mu * g3yi(:,4)
        qdi(:,12) =     rlm * u3 * g1yi(:,2) + rmu * u1 * g1yi(:,4)
     &            +    rlm * u3 * g2yi(:,3) + rmu * u2 * g2yi(:,4)
     &            +    rmu * u1 * g3yi(:,2) + rmu * u2 * g3yi(:,3)
     &            + rlm2mu * u3 * g3yi(:,4)
     &            +    con      * g3yi(:,5)

c
c
c.... assemble contribution of qdi to ql,i.e., contribution to
c     each element node
c
        do i=1,nshl
           ql(:,i,1 ) = ql(:,i,1 )+ shape(:,i)*WdetJ*qdi(:,1 )
           ql(:,i,2 ) = ql(:,i,2 )+ shape(:,i)*WdetJ*qdi(:,2 )
           ql(:,i,3 ) = ql(:,i,3 )+ shape(:,i)*WdetJ*qdi(:,3 )
           ql(:,i,4 ) = ql(:,i,4 )+ shape(:,i)*WdetJ*qdi(:,4 )
           ql(:,i,5 ) = ql(:,i,5 )+ shape(:,i)*WdetJ*qdi(:,5 )
           ql(:,i,6 ) = ql(:,i,6 )+ shape(:,i)*WdetJ*qdi(:,6 )
           ql(:,i,7 ) = ql(:,i,7 )+ shape(:,i)*WdetJ*qdi(:,7 )
           ql(:,i,8 ) = ql(:,i,8 )+ shape(:,i)*WdetJ*qdi(:,8 )
           ql(:,i,9 ) = ql(:,i,9 )+ shape(:,i)*WdetJ*qdi(:,9 )
           ql(:,i,10) = ql(:,i,10)+ shape(:,i)*WdetJ*qdi(:,10)
           ql(:,i,11) = ql(:,i,11)+ shape(:,i)*WdetJ*qdi(:,11)
           ql(:,i,12) = ql(:,i,12)+ shape(:,i)*WdetJ*qdi(:,12)
        enddo
c
c.... compute and assemble the element contribution to the lumped
c     mass matrix
c
c
c.... row sum technique
c
        if ( idiff == 1 ) then
           do i=1,nshl
              rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ
           enddo
        endif
c
c.... "special lumping technique" (Hughes p. 445)
c
        if ( idiff == 3 ) then
           shpsum = zero
           do i=1,nshl
              shpsum = shpsum + shape(:,i)*shape(:,i)
              rmassl(:,i)=rmassl(:,i)+shape(:,i)*shape(:,i)*WdetJ
           enddo
           alph1 = alph1+WdetJ
           alph2 = alph2+shpsum*WdetJ
        endif
      endif  ! end of idiff=1 .or. 3
c
      if(isurf .eq. 1) then
c
c.... initialize
c
        g1yti   = zero
        g2yti   = zero
        g3yti   = zero
c
c.... calculate the integration variables necessary for the
c     formation of q
c
c.... compute the global gradient of Yt-variables, assuming 6th entry as
c.... the phase indicator function
c
c  Yt_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(int) Yta)
c
        do n = 1, nshl
          g1yti(:)  = g1yti(:)  + shg(:,n,1) * ycl(:,n,6)
          g2yti(:)  = g2yti(:)  + shg(:,n,2) * ycl(:,n,6)
          g3yti(:)  = g3yti(:)  + shg(:,n,3) * ycl(:,n,6)
        enddo
c
c    computing N_{b}*N_{a,x_i)*yta*WdetJ
c
        do i=1,nshl
           ql(:,i,idflow+1)  = ql(:,i,idflow+1)
     &                       + shape(:,i)*WdetJ*g1yti
           ql(:,i,idflow+2)  = ql(:,i,idflow+2)
     &                       + shape(:,i)*WdetJ*g2yti
           ql(:,i,idflow+3)  = ql(:,i,idflow+3)
     &                       + shape(:,i)*WdetJ*g3yti
           rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ
        enddo
      endif  !end of the isurf
c
c.... end of the loop over integration points
c
      enddo
c
c.... normalize the mass matrix for idiff == 3
c
      if ( idiff == 3 ) then
         do i=1,nshl
            rmassl(:,i) = rmassl(:,i)*alph1/alph2
         enddo
      endif

      ttim(33) = ttim(33) + secs(0.0)

c
c.... return
c
       return
       end