xref: /phasta/phSolver/compressible/e3dc.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)

	subroutine e3DC (g1yi,   g2yi,   g3yi,   A0,     raLS,
     &			 rtLS,   giju,   DC,     ri,
     &                   rmi,    stiff, A0DC)
c
c----------------------------------------------------------------------
c
c This routine calculates the contribution of the Discontinuity-
c Capturing operator to RHS and preconditioner.
c
c  g1yi   (nflow,npro)           : grad-y in direction 1
c  g2yi   (nflow,npro)           : grad-y in direction 2
c  g3yi   (nflow,npro)           : grad-y in direction 3
c  A0     (nsymdf,npro)          : A0 matrix (Symm. storage)
c  raLS   (npro)                 : square of LS residual (A0inv norm)
c  rtLS   (npro)                 : square of LS residual (Tau norm)
c  giju    (6,npro)              : metric matrix
c  DC     (ngauss,npro)          : discontinuity-capturing factor
c  intp				 : integration point number
c
c output:
c  ri     (nflow*(nsd+1),npro)   : partial residual
c  rmi    (nflow*(nsd+1),npro)   : partial modified residual
c  stiff  (nsymdf,6,npro)       : diffusivity matrix
c  DC     (npro)                : discontinuity-capturing factor
c
c
c Zdenek Johan, Summer 1990. (Modified from e2dc.f)
c Zdenek Johan, Winter 1991. (Recoded)
c Zdenek Johan, Winter 1991. (Fortran 90)
c----------------------------------------------------------------------
c
	include "common.h"
c
        dimension g1yi(npro,nflow),          g2yi(npro,nflow),
     &            g3yi(npro,nflow),          A0(npro,5,5),
     &            raLS(npro),                rtLS(npro),
     &            giju(npro,6),              DC(npro,ngauss),
     &            ri(npro,nflow*(nsd+1)),    rmi(npro,nflow*(nsd+1)),
     &            stiff(npro,3*nflow,3*nflow),dtmp(npro)
c

        dimension ggyi(npro,nflow),         gAgyi(npro,15),
     &            gnorm(npro),              A0gyi(npro,15),
     &            yiA0DCyj(npro,6),         A0DC(npro,4)
c
c ... -----------------------> initialize <----------------------------
c
        A0gyi    = zero
        gAgyi    = zero
        yiA0DCyj = zero
        DC       = zero
c.... ----------------------->  global gradient  <----------------------
c
c.... calculate (A0 y_,j) --> A0gyi
c
c  A0 Y_{,1}
c
        A0gyi( :,1) = A0(:,1,1)*g1yi(:,1)
     &              + A0(:,1,2)*g1yi(:,2)
     &              + A0(:,1,3)*g1yi(:,3)
     &              + A0(:,1,4)*g1yi(:,4)
     &              + A0(:,1,5)*g1yi(:,5)
        A0gyi( :,2) = A0(:,2,1)*g1yi(:,1)
     &              + A0(:,2,2)*g1yi(:,2)
     &              + A0(:,2,3)*g1yi(:,3)
     &              + A0(:,2,4)*g1yi(:,4)
     &              + A0(:,2,5)*g1yi(:,5)
        A0gyi( :,3) = A0(:,3,1)*g1yi(:,1)
     &              + A0(:,3,2)*g1yi(:,2)
     &              + A0(:,3,3)*g1yi(:,3)
     &              + A0(:,3,4)*g1yi(:,4)
     &              + A0(:,3,5)*g1yi(:,5)
        A0gyi( :,4) = A0(:,4,1)*g1yi(:,1)
     &              + A0(:,4,2)*g1yi(:,2)
     &              + A0(:,4,3)*g1yi(:,3)
     &              + A0(:,4,4)*g1yi(:,4)
     &              + A0(:,4,5)*g1yi(:,5)
        A0gyi( :,5) = A0(:,5,1)*g1yi(:,1)
     &              + A0(:,5,2)*g1yi(:,2)
     &              + A0(:,5,3)*g1yi(:,3)
     &              + A0(:,5,4)*g1yi(:,4)
     &              + A0(:,5,5)*g1yi(:,5)
c
c  A0 Y_{,2}
c
        A0gyi( :,6) = A0(:,1,1)*g2yi(:,1)
     &              + A0(:,1,2)*g2yi(:,2)
     &              + A0(:,1,3)*g2yi(:,3)
     &              + A0(:,1,4)*g2yi(:,4)
     &              + A0(:,1,5)*g2yi(:,5)
        A0gyi( :,7) = A0(:,2,1)*g2yi(:,1)
     &              + A0(:,2,2)*g2yi(:,2)
     &              + A0(:,2,3)*g2yi(:,3)
     &              + A0(:,2,4)*g2yi(:,4)
     &              + A0(:,2,5)*g2yi(:,5)
        A0gyi( :,8) = A0(:,3,1)*g2yi(:,1)
     &              + A0(:,3,2)*g2yi(:,2)
     &              + A0(:,3,3)*g2yi(:,3)
     &              + A0(:,3,4)*g2yi(:,4)
     &              + A0(:,3,5)*g2yi(:,5)
        A0gyi( :,9) = A0(:,4,1)*g2yi(:,1)
     &              + A0(:,4,2)*g2yi(:,2)
     &              + A0(:,4,3)*g2yi(:,3)
     &              + A0(:,4,4)*g2yi(:,4)
     &              + A0(:,4,5)*g2yi(:,5)
        A0gyi(:,10) = A0(:,5,1)*g2yi(:,1)
     &              + A0(:,5,2)*g2yi(:,2)
     &              + A0(:,5,3)*g2yi(:,3)
     &              + A0(:,5,4)*g2yi(:,4)
     &              + A0(:,5,5)*g2yi(:,5)
c
c  A0 Y_{,3}
c
        A0gyi(:,11) = A0(:,1,1)*g3yi(:,1)
     &              + A0(:,1,2)*g3yi(:,2)
     &              + A0(:,1,3)*g3yi(:,3)
     &              + A0(:,1,4)*g3yi(:,4)
     &              + A0(:,1,5)*g3yi(:,5)
        A0gyi(:,12) = A0(:,2,1)*g3yi(:,1)
     &              + A0(:,2,2)*g3yi(:,2)
     &              + A0(:,2,3)*g3yi(:,3)
     &              + A0(:,2,4)*g3yi(:,4)
     &              + A0(:,2,5)*g3yi(:,5)
        A0gyi(:,13) = A0(:,3,1)*g3yi(:,1)
     &              + A0(:,3,2)*g3yi(:,2)
     &              + A0(:,3,3)*g3yi(:,3)
     &              + A0(:,3,4)*g3yi(:,4)
     &              + A0(:,3,5)*g3yi(:,5)
        A0gyi(:,14) = A0(:,4,1)*g3yi(:,1)
     &              + A0(:,4,2)*g3yi(:,2)
     &              + A0(:,4,3)*g3yi(:,3)
     &              + A0(:,4,4)*g3yi(:,4)
     &              + A0(:,4,5)*g3yi(:,5)
        A0gyi(:,15) = A0(:,5,1)*g3yi(:,1)
     &              + A0(:,5,2)*g3yi(:,2)
     &              + A0(:,5,3)*g3yi(:,3)
     &              + A0(:,5,4)*g3yi(:,4)
     &              + A0(:,5,5)*g3yi(:,5)
c
c.... calculate (giju A0 y_,j) --> gAgyi
c

        gAgyi( :,1) = giju(:,1)*A0gyi( :,1)
     &              + giju(:,4)*A0gyi( :,6)
     &              + giju(:,5)*A0gyi(:,11)

        gAgyi( :,2) = giju(:,1)*A0gyi( :,2)
     &              + giju(:,4)*A0gyi( :,7)
     &              + giju(:,5)*A0gyi(:,12)

	gAgyi( :,3) = giju(:,1)*A0gyi( :,3)
     &              + giju(:,4)*A0gyi( :,8)
     &              + giju(:,5)*A0gyi(:,13)

	gAgyi( :,4) = giju(:,1)*A0gyi( :,4)
     &              + giju(:,4)*A0gyi( :,9)
     &              + giju(:,5)*A0gyi(:,14)

	gAgyi( :,5) = giju(:,1)*A0gyi( :,5)
     &              + giju(:,4)*A0gyi(:,10)
     &              + giju(:,5)*A0gyi(:,15)

	gAgyi( :,6) = giju(:,4)*A0gyi( :,1)
     &              + giju(:,2)*A0gyi( :,6)
     &              + giju(:,6)*A0gyi(:,11)

	gAgyi( :,7) = giju(:,4)*A0gyi( :,2)
     &              + giju(:,2)*A0gyi( :,7)
     &              + giju(:,6)*A0gyi(:,12)

	gAgyi( :,8) = giju(:,4)*A0gyi( :,3)
     &              + giju(:,2)*A0gyi( :,8)
     &              + giju(:,6)*A0gyi(:,13)

	gAgyi( :,9) = giju(:,4)*A0gyi( :,4)
     &              + giju(:,2)*A0gyi( :,9)
     &              + giju(:,6)*A0gyi(:,14)

	gAgyi(:,10) = giju(:,4)*A0gyi( :,5)
     &              + giju(:,2)*A0gyi(:,10)
     &              + giju(:,6)*A0gyi(:,15)

	gAgyi(:,11) = giju(:,5)*A0gyi( :,1)
     &              + giju(:,6)*A0gyi( :,6)
     &              + giju(:,3)*A0gyi(:,11)

	gAgyi(:,12) = giju(:,5)*A0gyi( :,2)
     &              + giju(:,6)*A0gyi( :,7)
     &              + giju(:,3)*A0gyi(:,12)

	gAgyi(:,13) = giju(:,5)*A0gyi( :,3)
     &              + giju(:,6)*A0gyi( :,8)
     &              + giju(:,3)*A0gyi(:,13)

	gAgyi(:,14) = giju(:,5)*A0gyi( :,4)
     &              + giju(:,6)*A0gyi( :,9)
     &              + giju(:,3)*A0gyi(:,14)

	gAgyi(:,15) = giju(:,5)*A0gyi( :,5)
     &              + giju(:,6)*A0gyi(:,10)
     &              + giju(:,3)*A0gyi(:,15)
c
c... the denominator term of the DC factor
c... evaluation of the term  Y,i.A0DC Y,j
c
        yiA0DCyj(:,1) = A0DC(:,1)*g1yi(:,1)**2
     &                + two*g1yi(:,1)*A0DC(:,2)*g1yi(:,5)
     &                + A0DC(:,3)*g1yi(:,2)**2
     &                + A0DC(:,3)*g1yi(:,3)**2
     &                + A0DC(:,3)*g1yi(:,4)**2
     &                + A0DC(:,4)*g1yi(:,5)**2

        yiA0DCyj(:,2) = A0DC(:,1)*g2yi(:,1)**2
     &                + two*g2yi(:,1)*A0DC(:,2)*g2yi(:,5)
     &                + A0DC(:,3)*g2yi(:,2)**2
     &                + A0DC(:,3)*g2yi(:,3)**2
     &                + A0DC(:,3)*g2yi(:,4)**2
     &                + A0DC(:,4)*g2yi(:,5)**2

        yiA0DCyj(:,3) = A0DC(:,1)*g3yi(:,1)**2
     &                + two*g3yi(:,1)*A0DC(:,2)*g3yi(:,5)
     &                + A0DC(:,3)*g3yi(:,2)**2
     &                + A0DC(:,3)*g3yi(:,3)**2
     &                + A0DC(:,3)*g3yi(:,4)**2
     &                + A0DC(:,4)*g3yi(:,5)**2

        yiA0DCyj(:,4) = g1yi(:,1)*A0DC(:,1)*g2yi(:,1)
     &                + g1yi(:,1)*A0DC(:,2)*g2yi(:,5)
     &                + g1yi(:,2)*A0DC(:,3)*g2yi(:,2)
     &                + g1yi(:,3)*A0DC(:,3)*g2yi(:,3)
     &                + g1yi(:,4)*A0DC(:,3)*g2yi(:,4)
     &                + g1yi(:,5)*A0DC(:,2)*g2yi(:,1)
     &                + g1yi(:,5)*A0DC(:,4)*g2yi(:,5)

        yiA0DCyj(:,5) = g1yi(:,1)*A0DC(:,1)*g3yi(:,1)
     &                + g1yi(:,1)*A0DC(:,2)*g3yi(:,5)
     &                + g1yi(:,2)*A0DC(:,3)*g3yi(:,2)
     &                + g1yi(:,3)*A0DC(:,3)*g3yi(:,3)
     &                + g1yi(:,4)*A0DC(:,3)*g3yi(:,4)
     &                + g1yi(:,5)*A0DC(:,2)*g3yi(:,1)
     &                + g1yi(:,5)*A0DC(:,4)*g3yi(:,5)

        yiA0DCyj(:,6) = g2yi(:,1)*A0DC(:,1)*g3yi(:,1)
     &                + g2yi(:,1)*A0DC(:,2)*g3yi(:,5)
     &                + g2yi(:,2)*A0DC(:,3)*g3yi(:,2)
     &                + g2yi(:,3)*A0DC(:,3)*g3yi(:,3)
     &                + g2yi(:,4)*A0DC(:,3)*g3yi(:,4)
     &                + g2yi(:,5)*A0DC(:,2)*g3yi(:,1)
     &                + g2yi(:,5)*A0DC(:,4)*g3yi(:,5)
c
c.... ------------------------->  DC factor  <--------------------------
c
c	if ((ires .ne. 2) .or. (Jactyp .eq. 1)) then
c
c.... calculate 2-norm of Grad-local-V with respect to A0
c
c.... DC-mallet
c
	  if (iDC .eq. 1) then
c
	    fact = one
	    if (ipord .eq. 2)  fact = 0.9
	    if (ipord .eq. 3) fact = 0.75

c
            gnorm = one / (
     &              giju(:,1)*yiA0DCyj(:,1)
     &            + two*giju(:,4)*yiA0DCyj(:,4)
     &            + two*giju(:,5)*yiA0DCyj(:,5)
     &            + giju(:,2)*yiA0DCyj(:,2)
     &            + two*giju(:,6)*yiA0DCyj(:,6)
     &            + giju(:,3)*yiA0DCyj(:,3)
     &            + epsM  )
c
c	    DC(:,intp)=dim((fact*sqrt(raLS*gnorm)),(rtLS*gnorm))
	    DC(:,intp)=max(zero,(fact*sqrt(raLS*gnorm))-(rtLS*gnorm))
c
c.... flop count
c
!	    flops = flops + 46*npro
c
	  endif
c
c.... DC-quadratic
c
	  if (iDC .eq. 2) then
c
            gnorm = one / (
     &              giju(:,1)*yiA0DCyj(:,1)
     &            + two*giju(:,4)*yiA0DCyj(:,4)
     &            + two*giju(:,5)*yiA0DCyj(:,5)
     &            + giju(:,2)*yiA0DCyj(:,2)
     &            + two*giju(:,6)*yiA0DCyj(:,6)
     &            + giju(:,3)*yiA0DCyj(:,3)
     &            + epsM  )

c
	    DC(:,intp) = two * rtLS * gnorm
c
c.... flop count
c
!	    flops = flops + 36*npro
c
	  endif
c
c.... DC-min
c
	  if (iDC .eq. 3) then
c
	    fact = one
	    if (ipord .eq. 2)  fact = pt5
c
            gnorm = one / (
     &              giju(:,1)*yiA0DCyj(:,1)
     &            + two*giju(:,4)*yiA0DCyj(:,4)
     &            + two*giju(:,5)*yiA0DCyj(:,5)
     &            + giju(:,2)*yiA0DCyj(:,2)
     &            + two*giju(:,6)*yiA0DCyj(:,6)
     &            + giju(:,3)*yiA0DCyj(:,3)
     &            + epsM  )

c
c	    DC(:,intp) = min( dim(fact * sqrt(raLS * gnorm),
	    DC(:,intp) = min( max(zero,fact * sqrt(raLS * gnorm)-
     &                       rtLS * gnorm), two * rtLS * gnorm )
c
c.... flop count
c
!	    flops = flops + 48*npro
c
	  endif
c
c	endif
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... add the contribution of DC to ri and/or rmi
c
c.... ires = 1 or 3
c
	if ((ires .eq. 1) .or. (ires .eq. 3)) then
c
	  ri ( :,1) = ri ( :,1) + DC(:,intp) * gAgyi( :,1)
	  rmi( :,1) = rmi( :,1) + DC(:,intp) * gAgyi( :,1)
	  ri ( :,2) = ri ( :,2) + DC(:,intp) * gAgyi( :,2)
	  rmi( :,2) = rmi( :,2) + DC(:,intp) * gAgyi( :,2)
	  ri ( :,3) = ri ( :,3) + DC(:,intp) * gAgyi( :,3)
	  rmi( :,3) = rmi( :,3) + DC(:,intp) * gAgyi( :,3)
	  ri ( :,4) = ri ( :,4) + DC(:,intp) * gAgyi( :,4)
	  rmi( :,4) = rmi( :,4) + DC(:,intp) * gAgyi( :,4)
	  ri ( :,5) = ri ( :,5) + DC(:,intp) * gAgyi( :,5)
	  rmi( :,5) = rmi( :,5) + DC(:,intp) * gAgyi( :,5)
c
	  ri ( :,6) = ri ( :,6) + DC(:,intp) * gAgyi( :,6)
	  rmi( :,6) = rmi( :,6) + DC(:,intp) * gAgyi( :,6)
	  ri ( :,7) = ri ( :,7) + DC(:,intp) * gAgyi( :,7)
	  rmi( :,7) = rmi( :,7) + DC(:,intp) * gAgyi( :,7)
	  ri ( :,8) = ri ( :,8) + DC(:,intp) * gAgyi( :,8)
	  rmi( :,8) = rmi( :,8) + DC(:,intp) * gAgyi( :,8)
	  ri ( :,9) = ri ( :,9) + DC(:,intp) * gAgyi( :,9)
	  rmi( :,9) = rmi( :,9) + DC(:,intp) * gAgyi( :,9)
	  ri (:,10) = ri (:,10) + DC(:,intp) * gAgyi(:,10)
	  rmi(:,10) = rmi(:,10) + DC(:,intp) * gAgyi(:,10)
c
	  ri (:,11) = ri (:,11) + DC(:,intp) * gAgyi(:,11)
	  rmi(:,11) = rmi(:,12) + DC(:,intp) * gAgyi(:,12)
	  ri (:,12) = ri (:,12) + DC(:,intp) * gAgyi(:,12)
	  rmi(:,12) = rmi(:,12) + DC(:,intp) * gAgyi(:,12)
	  ri (:,13) = ri (:,13) + DC(:,intp) * gAgyi(:,13)
	  rmi(:,13) = rmi(:,13) + DC(:,intp) * gAgyi(:,13)
	  ri (:,14) = ri (:,14) + DC(:,intp) * gAgyi(:,14)
	  rmi(:,14) = rmi(:,14) + DC(:,intp) * gAgyi(:,14)
	  ri (:,15) = ri (:,15) + DC(:,intp) * gAgyi(:,15)
	  rmi(:,15) = rmi(:,15) + DC(:,intp) * gAgyi(:,15)
c
!	  flops = flops + 45*npro
c
	endif
c
c.... ires = 2
c
	if (ires .eq. 2) then
c
	  rmi( :,1) = rmi( :,1) + DC(:,intp) * gAgyi( :,1)
	  rmi( :,2) = rmi( :,2) + DC(:,intp) * gAgyi( :,2)
	  rmi( :,3) = rmi( :,3) + DC(:,intp) * gAgyi( :,3)
	  rmi( :,4) = rmi( :,4) + DC(:,intp) * gAgyi( :,4)
	  rmi( :,5) = rmi( :,5) + DC(:,intp) * gAgyi( :,5)
c
	  rmi( :,6) = rmi( :,6) + DC(:,intp) * gAgyi( :,6)
	  rmi( :,7) = rmi( :,7) + DC(:,intp) * gAgyi( :,7)
	  rmi( :,8) = rmi( :,8) + DC(:,intp) * gAgyi( :,8)
	  rmi( :,9) = rmi( :,9) + DC(:,intp) * gAgyi( :,9)
	  rmi(:,10) = rmi(:,10) + DC(:,intp) * gAgyi(:,10)
c
	  rmi(:,11) = rmi(:,11) + DC(:,intp) * gAgyi(:,11)
	  rmi(:,12) = rmi(:,12) + DC(:,intp) * gAgyi(:,12)
	  rmi(:,13) = rmi(:,13) + DC(:,intp) * gAgyi(:,13)
	  rmi(:,14) = rmi(:,14) + DC(:,intp) * gAgyi(:,14)
	  rmi(:,15) = rmi(:,15) + DC(:,intp) * gAgyi(:,15)
c
!	  flops = flops + 30*npro
c
	endif
c
c.... ------------------------->  Stiffness  <--------------------------
c
c.... add the contribution of DC to stiff
c
	if (iprec .eq. 1) then ! leave out of LHS, when called from itrres
	     nflow2=two*nflow
       do j = 1, nflow
          do i = 1, nflow
             dtmp(:)=A0(:,i,j)*DC(:,intp)
c
c.... add (DC g^1 A0) to stiff [1,1]
c
             stiff(:,i,j) = stiff(:,i,j)
     &                    + dtmp(:)*giju(:,1)
c
c.... add (DC g^1 A0) to stiff [1,2]
c

             stiff(:,i,j+nflow) = stiff(:,i,j+nflow)
     &                    + dtmp(:)*giju(:,4)
c
c.... add (DC g^1 A0) to stiff [1,3]
c

             stiff(:,i,j+nflow2) = stiff(:,i,j+nflow2)
     &                    + dtmp(:)*giju(:,5)

c.... add (DC g^1 A0) to stiff [2,1] (similarly below)
c

             stiff(:,i+nflow,j) = stiff(:,i+nflow,j)
     &                    + dtmp(:)*giju(:,4)

             stiff(:,i+nflow,j+nflow) = stiff(:,i+nflow,j+nflow)
     &                    + dtmp(:)*giju(:,2)

             stiff(:,i+nflow,j+nflow2) = stiff(:,i+nflow,j+nflow2)
     &                    + dtmp(:)*giju(:,6)

             stiff(:,i+nflow2,j) = stiff(:,i+nflow2,j)
     &                    + dtmp(:)*giju(:,5)

             stiff(:,i+nflow2,j+nflow) = stiff(:,i+nflow2,j+nflow)
     &                    + dtmp(:)*giju(:,6)

             stiff(:,i+nflow2,j+nflow2) = stiff(:,i+nflow2,j+nflow2)
     &                    + dtmp(:)*giju(:,3)
          enddo
       enddo
c
c.... flop count
c
!	  flops = flops + 210*npro
c
c.... end of stiffness
c
	endif
c
c.... return
c
	return
	end
c
        subroutine e3dcSclr ( g1yti,    g2yti,       g3yti,
     &                        A0t,      raLSt,       rTLSt,
     &                        DCt,      giju,
     &                        rti,      rmti,        stifft)
c
c
c----------------------------------------------------------------------
c
c This routine calculates the contribution of the Discontinuity-
c Capturing operator to RHS and preconditioner for the scalar solve.
c
c  g1yti   (nflow,npro)           : grad-y in direction 1
c  g2yti   (nflow,npro)           : grad-y in direction 2
c  g3yti   (nflow,npro)           : grad-y in direction 3
c  A0     (nsymdf,npro)          : A0 matrix (Symm. storage)
c  raLS   (npro)                 : square of LS residual (A0inv norm)
c  rtLS   (npro)                 : square of LS residual (Tau norm)
c  giju    (6,npro)              : metric matrix
c  DC     (ngauss,npro)          : discontinuity-capturing factor
c  intp				 : integration point number
c
c output:
c  ri     (nflow*(nsd+1),npro)   : partial residual
c  rmi    (nflow*(nsd+1),npro)   : partial modified residual
c  stiff  (nsymdf,6,npro)       : diffusivity matrix
c  DC     (npro)                : discontinuity-capturing factor
c
c
c Zdenek Johan, Summer 1990. (Modified from e2dc.f)
c Zdenek Johan, Winter 1991. (Recoded)
c Zdenek Johan, Winter 1991. (Fortran 90)
c----------------------------------------------------------------------
c
	include "common.h"
c
        dimension g1yti(npro),                g2yti(npro),
     &            g3yti(npro),                A0t(npro),
     &            raLSt(npro),                rtLSt(npro),
     &            giju(npro,6),               DCt(npro,ngauss),
     &            rti(npro,nsd+1),            rmti(npro,nsd+1),
     &            stifft(npro,nsd,nsd),       dtmp(npro)
c

        dimension ggyit(npro,nflow),        gAgyit(npro,3),
     &            gnormt(npro),             A0gyit(npro,3),
     &            yiA0DCyjt(npro,6),        A0DCt(npro)
c
c ... -----------------------> initialize <----------------------------
c
        A0gyit    = zero
        gAgyit    = zero
        yiA0DCyjt = zero
        DCt       = zero
        A0DCt    = A0t
c.... ----------------------->  global gradient  <----------------------
c
c.... calculate (A0 y_,j) --> A0gyit
c
c  A0 Y_{,1}
c
        A0gyit( :,1) = A0t(:)*g1yti(:)
c  A0 Y_{,2}
        A0gyit( :,2) = A0t(:)*g2yti(:)
c  A0 Y_{,3}
        A0gyit( :,3) = A0t(:)*g3yti(:)
c
c.... calculate (giju A0 y_,j) --> gAgyit
c

        gAgyit( :,1) = giju(:,1)*A0gyit( :,1)
     &               + giju(:,4)*A0gyit( :,2)
     &               + giju(:,5)*A0gyit( :,3)

        gAgyit( :,2) = giju(:,4)*A0gyit( :,1)
     &               + giju(:,2)*A0gyit( :,2)
     &               + giju(:,6)*A0gyit( :,3)

	gAgyit( :,3) = giju(:,5)*A0gyit( :,1)
     &               + giju(:,6)*A0gyit( :,2)
     &               + giju(:,3)*A0gyit( :,3)
c
c... the denominator term of the DC factor
c... evaluation of the term  Y,i.A0DC Y,j
c
        yiA0DCyjt(:,1) = A0DCt(:)*g1yti(:)**2
c
        yiA0DCyjt(:,2) = A0DCt(:)*g2yti(:)**2
c
        yiA0DCyjt(:,3) = A0DCt(:)*g3yti(:)**2
c
        yiA0DCyjt(:,4) = A0DCt(:)*g1yti(:)*g2yti(:)
c
        yiA0DCyjt(:,5) = A0DCt(:)*g1yti(:)*g3yti(:)
c
        yiA0DCyjt(:,6) = A0DCt(:)*g2yti(:)*g3yti(:)
c
c
c.... ------------------------->  DC factor  <--------------------------
c
c	if ((ires .ne. 2) .or. (Jactyp .eq. 1)) then
c
c.... calculate 2-norm of Grad-local-V with respect to A0
c
c.... DC-mallet
c
	  if (iDCsclr(1) .eq. 1) then
c
	    fact = one
	    if (ipord .eq. 2)  fact = 0.9
	    if (ipord .eq. 3) fact = 0.75

c
            gnormt = one / (
     &              giju(:,1)*yiA0DCyjt(:,1)
     &            + two*giju(:,4)*yiA0DCyjt(:,4)
     &            + two*giju(:,5)*yiA0DCyjt(:,5)
     &            + giju(:,2)*yiA0DCyjt(:,2)
     &            + two*giju(:,6)*yiA0DCyjt(:,6)
     &            + giju(:,3)*yiA0DCyjt(:,3)
     &            + epsM  )
c
c	    DCt(:,intp)=dim((fact*sqrt(raLSt*gnormt)),(rtLSt*gnormt))
	    DCt(:,intp)=max(zero,(fact*sqrt(raLSt*gnormt))-(rtLSt*gnormt))
c
c.... flop count
c
!	    flops = flops + 46*npro
c
	  endif
c
c.... DC-quadratic
c
	  if (iDCSclr(1) .eq. 2) then
c
            gnormt = one / (
     &              giju(:,1)*yiA0DCyjt(:,1)
     &            + two*giju(:,4)*yiA0DCyjt(:,4)
     &            + two*giju(:,5)*yiA0DCyjt(:,5)
     &            + giju(:,2)*yiA0DCyjt(:,2)
     &            + two*giju(:,6)*yiA0DCyjt(:,6)
     &            + giju(:,3)*yiA0DCyjt(:,3)
     &            + epsM  )

c
	    DCt(:,intp) = two * rtLSt * gnormt
c
c.... flop count
c
!	    flops = flops + 36*npro
c
	  endif
c
c.... DC-min
c
	  if (iDCSclr(1) .eq. 3) then
c
	    fact = one
	    if (ipord .eq. 2)  fact = pt5
c
            gnormt = one / (
     &              giju(:,1)*yiA0DCyjt(:,1)
     &            + two*giju(:,4)*yiA0DCyjt(:,4)
     &            + two*giju(:,5)*yiA0DCyjt(:,5)
     &            + giju(:,2)*yiA0DCyjt(:,2)
     &            + two*giju(:,6)*yiA0DCyjt(:,6)
     &            + giju(:,3)*yiA0DCyjt(:,3)
     &            + epsM  )

c
c	    DCt(:,intp) = min( dim(fact * sqrt(raLSt * gnormt),
	    DCt(:,intp) = min( max(zero,fact * sqrt(raLSt * gnormt)-
     &                       rtLSt * gnormt), two * rtLSt * gnormt )
c
c.... flop count
c
!	    flops = flops + 48*npro
c
	  endif
c
c	endif
c	DCt=DCt*two
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... add the contribution of DC to ri and/or rmi
c
c.... ires = 1 or 3
c
	if ((ires .eq. 1) .or. (ires .eq. 3)) then
c
	  rti ( :,1) = rti ( :,1) + DCt(:,intp) * gAgyit( :,1)
	  rmti( :,1) = rmti( :,1) + DCt(:,intp) * gAgyit( :,1)
	  rti ( :,2) = rti ( :,2) + DCt(:,intp) * gAgyit( :,2)
	  rmti( :,2) = rmti( :,2) + DCt(:,intp) * gAgyit( :,2)
	  rti ( :,3) = rti ( :,3) + DCt(:,intp) * gAgyit( :,3)
	  rmti( :,3) = rmti( :,3) + DCt(:,intp) * gAgyit( :,3)

c
!	  flops = flops + 45*npro
c
	endif
c
c.... ires = 2
c
	if (ires .eq. 2) then
c
	  rmti( :,1) = rmti( :,1) + DCt(:,intp) * gAgyit( :,1)
	  rmti( :,2) = rmti( :,2) + DCt(:,intp) * gAgyit( :,2)
	  rmti( :,3) = rmti( :,3) + DCt(:,intp) * gAgyit( :,3)

c
!	  flops = flops + 30*npro
c
	endif
c
c.... ------------------------->  Stiffness  <--------------------------
c
c.... add the contribution of DC to stiff
c$$$c
c	if (iprec .eq. 1) then !leave out of LHS, if called from itrres
	                       !anyway matrix free not implemented for scalar
             dtmp(:)=A0t(:)*DCt(:,intp)
c
c.... add (DC g^1 A0) to stifft [1,1]
c
             stifft(:,1,1) = stifft(:,1,1)
     &                    + dtmp(:)*giju(:,1)
c
c.... add (DC g^1 A0) to stifft [1,2]
c
             stifft(:,1,2) = stifft(:,1,2)
     &                    + dtmp(:)*giju(:,4)
c
c.... add (DC g^1 A0) to stifft [1,3]
c
             stifft(:,1,3) = stifft(:,1,3)
     &                    + dtmp(:)*giju(:,5)

c.... add (DC g^1 A0) to stifft [2,1] (similarly below)
c
             stifft(:,2,1) = stifft(:,2,1)
     &                    + dtmp(:)*giju(:,4)

             stifft(:,2,2) = stifft(:,2,2)
     &                    + dtmp(:)*giju(:,2)

             stifft(:,2,3) = stifft(:,2,3)
     &                    + dtmp(:)*giju(:,6)

             stifft(:,3,1) = stifft(:,3,1)
     &                    + dtmp(:)*giju(:,5)

             stifft(:,3,2) = stifft(:,3,2)
     &                    + dtmp(:)*giju(:,6)

             stifft(:,3,3) = stifft(:,3,3)
     &                    + dtmp(:)*giju(:,3)
c
c.... flop count
c
!	  flops = flops + 210*npro
c
c.... end of stiffness
c
c	endif
c
c.... return
c
	return
	end
c