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

        subroutine e3wmlt (shp,    shg,    WdetJ,
     &                     ri,     rmi,    rl,
     &                     rml,    stiff,  EGmass )
c
c----------------------------------------------------------------------
c
c Up to now most of the terms have not been multiplied by the
c shape function from the weight space.  Now that we have collected
c all the terms that the shape function (and its derivatives for the
c terms that were integrated by parts) multiplies, in this routine
c we carry out that multiplication. After these operations we will
c have this quadrature points contribution to the integral for the
c residual (i.e.  g^e_b(xi^l) D(xi^l) QW(xi^l)  where e is the element,
c b is the local node number on the element, l is the quadrature point
c D is the determinate of the Jacobian of the mapping, and QW is this
c quadrature points weight....WHEW).  When we add up all of the
c integration points we get G^e_b which we will assemble in the
c subroutine LOCAL to form G_B.
c
c This routine also forms the LHS contribution from the LS term
c and the diffusion term which has been partially constructed and
c placed in stiff.  Those familiar with elasticity might recognize
c this naming convention since this is like a stiffness matrix that
c if you had a linear problem would be calculated once and saved for
c all time.
c
c    ri:  LS, Diffusion, Convective and mass,
c   rmi:  LS, Diffusion and mass,
c stiff:  LS, Diffusion.
c
c input:
c  shp    (nshl)                 : element shape-functions
c  shg    (npro,nshl,nsd)        : element grad-shape-functions
c  WdetJ  (npro)                   : weighted Jacobian
c  ri     (npro,nflow*(nsd+1))      : partial residual
c  rmi    (npro,nflow*(nsd+1))      : partial modified residual
c  stiff  (npro,nsd*nflow,nsd*nflow) : stiffness matrix
c                                    ( K_ij + A_i tau A_j )
c
c output:
c  rl     (npro,nshl,nflow)       : residual
c  rml    (npro,nshl,nflow)       : modified residual
c  EGmass (npro,nedof,nedof)       : element LHS tangent mass matrix
c
c
c Zdenek Johan, Summer 1990.  (Modified from e2assm.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c Kenneth Jansen, Winter 1997  Prim. variables
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension shp(npro,nshl),
     &            shg(npro,nshl,nsd),
     &            WdetJ(npro),             ri(npro,nflow*(nsd+1)),
     &            rmi(npro,nflow*(nsd+1)), stiff(npro,3*nflow,3*nflow),
     &            rl(npro,nshl,nflow),     rml(npro,nshl,nflow),
     &            EGmass(npro,nedof,nedof)
c
        dimension shg1(npro),              shg2(npro),
     &            shg3(npro),              stif1(npro,nflow,nflow),
     &            stif2(npro,nflow,nflow), stif3(npro,nflow,nflow)

        ttim(29) = ttim(29) - secs(0.0)
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... add spatial contribution to rl and rml
c
c.... ires = 1 or 3
c
        if ((ires .eq. 1) .or. (ires .eq. 3)) then
c
          do i = 1, nshl
            rl(:,i,1) = rl(:,i,1) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 1) + shg(:,i,2) * ri(:, 6)
     &                                        + shg(:,i,3) * ri(:,11) )
            rl(:,i,2) = rl(:,i,2) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 2) + shg(:,i,2) * ri(:, 7)
     &                                        + shg(:,i,3) * ri(:,12) )
            rl(:,i,3) = rl(:,i,3) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 3) + shg(:,i,2) * ri(:, 8)
     &                                        + shg(:,i,3) * ri(:,13) )
            rl(:,i,4) = rl(:,i,4) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 4) + shg(:,i,2) * ri(:, 9)
     &                                        + shg(:,i,3) * ri(:,14) )
            rl(:,i,5) = rl(:,i,5) + WdetJ * (
     &                  shg(:,i,1) * ri(:, 5) + shg(:,i,2) * ri(:,10)
     &                                        + shg(:,i,3) * ri(:,15) )
          enddo
c
          flops = flops + 36*nshl*npro
        endif
c
c.... ires = 2 or 3
c
        if ((ires .eq. 2) .or. (ires .eq. 3)) then
          do i = 1, nshl
            rml(:,i,1) = rml(:,i,1) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 1) + shg(:,i,2) * rmi(:, 6)
     &                 + shg(:,i,3) * rmi(:,11)
     &                   + shp(:,i) * rmi(:,16)  )
            rml(:,i,2) = rml(:,i,2) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 2) + shg(:,i,2) * rmi(:, 7)
     &                 + shg(:,i,3) * rmi(:,12)
     &                   + shp(:,i) * rmi(:,17)    )
            rml(:,i,3) = rml(:,i,3) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 3) + shg(:,i,2) * rmi(:, 8)
     &                 + shg(:,i,3) * rmi(:,13)
     &                   + shp(:,i) * rmi(:,18)    )
            rml(:,i,4) = rml(:,i,4) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 4) + shg(:,i,2) * rmi(:, 9)
     &                 + shg(:,i,3) * rmi(:,14)
     &                   + shp(:,i) * rmi(:,19)    )
            rml(:,i,5) = rml(:,i,5) + WdetJ * (
     &                   shg(:,i,1) * rmi(:, 5) + shg(:,i,2) * rmi(:,10)
     &                 + shg(:,i,3) * rmi(:,15)
     &                   + shp(:,i) * rmi(:,20)    )
          enddo
c
          flops = flops + (15+45*nshl)*npro
        endif
c
c.... add temporal contribution to rl
c
        if (ngauss .eq. 1 .and. nshl.eq.4) then !already ex. integ mass
         if(matflg(5,1).ge.1) then
          do i = 1, nshl
            rl(:,i,2) = rl(:,i,2) + shp(:,i) * WdetJ * ri(:,17)
            rl(:,i,3) = rl(:,i,3) + shp(:,i) * WdetJ * ri(:,18)
            rl(:,i,4) = rl(:,i,4) + shp(:,i) * WdetJ * ri(:,19)
            rl(:,i,5) = rl(:,i,5) + shp(:,i) * WdetJ * ri(:,20)
          enddo
         endif
       else   !check for a body force
          do i = 1, nshl
            rl(:,i,1) = rl(:,i,1) + shp(:,i) * WdetJ * ri(:,16)
            rl(:,i,2) = rl(:,i,2) + shp(:,i) * WdetJ * ri(:,17)
            rl(:,i,3) = rl(:,i,3) + shp(:,i) * WdetJ * ri(:,18)
            rl(:,i,4) = rl(:,i,4) + shp(:,i) * WdetJ * ri(:,19)
            rl(:,i,5) = rl(:,i,5) + shp(:,i) * WdetJ * ri(:,20)
          enddo

          flops = flops + 11*nshl*npro
       endif

c
c.... ---------------------------->  LHS  <----------------------------
c
       if (lhs .eq. 1) then
c
c.... loop through columns (nodes j)
c
          do j = 1, nshl
             j0 = nflow * (j - 1)
c
c.... set up factors
c
             shg1 = WdetJ * shg(:,j,1)
             shg2 = WdetJ * shg(:,j,2)
             shg3 = WdetJ * shg(:,j,3)
c
c.... loop through d.o.f.'s
c
             do jdof = 1, nflow
                do idof = 1, nflow
                   idof2 = idof  + nflow
                   jdof2 = jdof  + nflow

                   idof3 = idof2 + nflow
                   jdof3 = jdof2 + nflow
c
c.... calculate weighted stiffness matrix (first part)
c
                   stif1(:,idof,jdof) = shg1 * stiff(:,idof,jdof)
     &                                + shg2 * stiff(:,idof,jdof2)
     &                                + shg3 * stiff(:,idof,jdof3)
                   stif2(:,idof,jdof) = shg1 * stiff(:,idof2,jdof)
     &                                + shg2 * stiff(:,idof2,jdof2)
     &                                + shg3 * stiff(:,idof2,jdof3)
                   stif3(:,idof,jdof) = shg1 * stiff(:,idof3,jdof)
     &                                + shg2 * stiff(:,idof3,jdof2)
     &                                + shg3 * stiff(:,idof3,jdof3)
                enddo
             enddo
c
c.... loop through rows (nodes i)
c
             do i = 1, nshl
                i0 = nflow * (i - 1)
c
c.... add contribution of stiffness to EGmass
c
                do jdof = 1, nflow
                   EGmass(:,i0+1,j0+jdof) = EGmass(:,i0+1,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,1,jdof)
     &                  + shg(:,i,2) * stif2(:,1,jdof)
     &                  + shg(:,i,3) * stif3(:,1,jdof)
                   EGmass(:,i0+2,j0+jdof) = EGmass(:,i0+2,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,2,jdof)
     &                  + shg(:,i,2) * stif2(:,2,jdof)
     &                  + shg(:,i,3) * stif3(:,2,jdof)
                   EGmass(:,i0+3,j0+jdof) = EGmass(:,i0+3,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,3,jdof)
     &                  + shg(:,i,2) * stif2(:,3,jdof)
     &                  + shg(:,i,3) * stif3(:,3,jdof)
                   EGmass(:,i0+4,j0+jdof) = EGmass(:,i0+4,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,4,jdof)
     &                  + shg(:,i,2) * stif2(:,4,jdof)
     &                  + shg(:,i,3) * stif3(:,4,jdof)
                   EGmass(:,i0+5,j0+jdof) = EGmass(:,i0+5,j0+jdof)
     &                  + shg(:,i,1) * stif1(:,5,jdof)
     &                  + shg(:,i,2) * stif2(:,5,jdof)
     &                  + shg(:,i,3) * stif3(:,5,jdof)
                enddo
c
c.... end loop on rows
c
             enddo
c
c.... end loop on columns
c
          enddo
c
c.... end left hand side assembly
c
       endif

       ttim(29) = ttim(29) + secs(0.0)
c
c.... return
c
        return
        end
c
c
c
        subroutine e3wmltSclr(shp,    shg,    WdetJ,
     &                        rti,    rtl,
     &                        stifft, EGmasst )
c
c----------------------------------------------------------------------
c
c Up to now most of the terms have not been multiplied by the
c shape function from the weight space.  Now that we have collected
c all the terms that the shape function (and its derivatives for the
c terms that were integrated by parts) multiplies, in this routine
c we carry out that multiplication. After these operations we will
c have this quadrature points contribution to the integral for the
c residual (i.e.  g^e_b(xi^l) D(xi^l) QW(xi^l)  where e is the element,
c b is the local node number on the element, l is the quadrature point
c D is the determinate of the Jacobian of the mapping, and QW is this
c quadrature points weight....WHEW).  When we add up all of the
c integration points we get G^e_b which we will assemble in the
c subroutine LOCAL to form G_B.
c
c This routine also forms the LHS contribution from the LS term
c and the diffusion term which has been partially constructed and
c placed in stiff.  Those familiar with elasticity might recognize
c this naming convention since this is like a stiffness matrix that
c if you had a linear problem would be calculated once and saved for
c all time.
c
c    ri:  LS, Diffusion, Convective and mass,
c stiff:  LS, Diffusion.
c
c input:
c  shp     (npro,nshl)            : element shape-functions
c  shg     (npro,nshl,nsd)        : element grad-shape-functions
c  WdetJ   (npro)                 : weighted Jacobian
c  rti     (npro,nsd+1)           : partial residual
c  stifft  (npro,nsd,nsd)         : stiffness matrix
c                                    ( K_ij + A_i tau A_j )
c
c output:
c  rtl     (npro,nshl)            : residual  (will end up being G^e_b)
c  EGmasst (npro,nshape,nshape)   : element LHS tangent mass matrix
c                                  (will end up being dG^e_b/dY_a)
c
c
c Zdenek Johan, Summer 1990.  (Modified from e2assm.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c Kenneth Jansen, Winter 1997  Prim. variables
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension shp(npro,nshl),            shg(npro,nshl,nsd),
     &            WdetJ(npro),               rti(npro,nsd+1),
     &            rmti(npro,nsd+1),          stifft(npro,nsd,nsd),
     &            rtl(npro,nshl),            rmtl(npro,nshl),
     &            EGmasst(npro,nshape,nshape)
c
        dimension shg1(npro),                shg2(npro),
     &            shg3(npro),                stift1(npro),
     &            stift2(npro),              stift3(npro)

        ttim(29) = ttim(29) - tmr(0.0)
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... add spatial contribution to rl and rml
c
c.... ires = 1 or 3
c
        if ((ires .eq. 1) .or. (ires .eq. 3)) then
c
          do i = 1, nshl
            rtl(:,i) = rtl(:,i) + WdetJ * (
     &                  shg(:,i,1) * rti(:,1) + shg(:,i,2) * rti(:,2)
     &                                        + shg(:,i,3) * rti(:,3) )

          enddo
          flops = flops + 36*nshl*npro
        endif
c
c.... ires = 2 or 3
c
c        if ((ires .eq. 2) .or. (ires .eq. 3)) then
c          do i = 1, nshl
c            rml(:,i,1) = rml(:,i,1) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 1) + shg(:,i,2) * rmi(:, 6)
c     &                 + shg(:,i,3) * rmi(:,11) + shp(:,i) * rmi(:,16) )
c            rml(:,i,2) = rml(:,i,2) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 2) + shg(:,i,2) * rmi(:, 7)
c     &                 + shg(:,i,3) * rmi(:,12) + shp(:,i) * rmi(:,17) )
c            rml(:,i,3) = rml(:,i,3) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 3) + shg(:,i,2) * rmi(:, 8)
c     &                 + shg(:,i,3) * rmi(:,13) + shp(:,i) * rmi(:,18) )
c            rml(:,i,4) = rml(:,i,4) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 4) + shg(:,i,2) * rmi(:, 9)
c     &                 + shg(:,i,3) * rmi(:,14) + shp(:,i) * rmi(:,19) )
c            rml(:,i,5) = rml(:,i,5) + WdetJ * (
c     &                   shg(:,i,1) * rmi(:, 5) + shg(:,i,2) * rmi(:,10)
c     &                 + shg(:,i,3) * rmi(:,15) + shp(:,i) * rmi(:,20) )
c          enddo
c
c          flops = flops + (15+45*nshl)*npro
c        endif

c
c.... add temporal contribution to rl
c

          do i = 1, nshl
            rtl(:,i) = rtl(:,i) + shp(:,i) * WdetJ * rti(:,4)
          enddo

          flops = flops + 11*nshl*npro

c
c.... ---------------------------->  LHS  <----------------------------
c
       if (lhs .eq. 1) then
c
c.... loop through columns (nodes j)
c
          do j = 1, nshl
c
c.... set up factors
c
             shg1 = WdetJ * shg(:,j,1)
             shg2 = WdetJ * shg(:,j,2)
             shg3 = WdetJ * shg(:,j,3)
c
c.... calculate weighted stiffness matrix (first part)
c
                   stift1(:) = shg1 * stifft(:,1,1)
     &                                + shg2 * stifft(:,1,2)
     &                                + shg3 * stifft(:,1,3)
                   stift2(:) = shg1 * stifft(:,2,1)
     &                                + shg2 * stifft(:,2,2)
     &                                + shg3 * stifft(:,2,3)
                   stift3(:) = shg1 * stifft(:,3,1)
     &                                + shg2 * stifft(:,3,2)
     &                                + shg3 * stifft(:,3,3)
c
c.... loop through rows (nodes i)
c
             do i = 1, nshl
c
c.... add contribution of stiffness to EGmass
c
                   EGmasst(:,i,j) = EGmasst(:,i,j)
     &                  + shg(:,i,1) * stift1(:)
     &                  + shg(:,i,2) * stift2(:)
     &                  + shg(:,i,3) * stift3(:)

c
c.... end loop on rows
c
             enddo
c
c.... end loop on columns
c
          enddo
c
c.... end left hand side assembly
c
       endif

       ttim(29) = ttim(29) + tmr(0.0)
c
c.... return
c
        return
        end