C tot_energy_2d_1.f

      subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,ni,layer,nndi,
     * nnshear,jpltcd)
C
C  This subroutine computes the area and total strain energy of a
C   4 noded 2D element mesh.
C       Post code       -1      Area
C       Post code       -2      Strain energy density
C       Post code       -3      Total strain energy
c* * * * * *
C
      implicit real*8 (a-h,o-z)
      dimension x(2,4),xintp(2,4),f(4),dfdx(4),dfdy(4),dfdz(4)
      dimension s(1),etot(1),eplas(1),ecreep(1),sp(1)
      include '../common/concom'
      include '../common/creeps'
      include '../common/elmcom'
      include '../common/matdat'
      include '../common/heat'
      include '../common/lass'
      include '../common/array2'
      include '../common/arrays'
      include '../common/space'

        include '../common/nzro1'
        include '../common/strnen'

	common/myvar/total,incold

      if (inc .gt. incold) then
       write(75,*)"Total strain energy during inc.",incold," is",total
       total=0.
       write(75,*)
      endif
      incold = inc

C **********************************
        lms = nnx-1
        lms2 = lms*neqst
        mms2 = lms2
        lofr = (m-1)*nelstr
        jtoten = itoten + lofr + mms2
        totde = vars(jtoten)

C **********************************

	la1=icoord
	call ppp003(f,dfdx,dfdy,dfdz,vars(la1),dvol,xintp,m,ni,ierr)
	if (jpltcd .eq. 1) v=4.*dvol
	if (jpltcd .eq. 2) v=totde
	if (jpltcd .eq. 3) then
	  v=4*dvol*totde
	  total = total+v
        write(76,*)"Increment",inc," Total Energy",total
        endif
	
	return
	end


      subroutine ppp003(f,dfdx,dfdy,dfdz,x,dvol,xintp,m,nn,ierr)
      implicit real*8 (a-h,o-z)
c
      dimension x(2,4),xintp(2,4),f(4),dfdx(4),dfdy(4),dfdz(4),
     $          gausx(4),gausy(4),dfksi(4),dfeta(4)
c
      data gausx/-.577350269d0, .577350269d0,-.577350269d0,
     $ .577350269d0/
      data gausy/-.577350269d0,-.577350269d0, .577350269d0,
     $ .577350269d0/
c
      xsi = gausx(nn)
      eta = gausy(nn)
      xm=1.d0-xsi
      xp=1.d0+xsi
      ym=1.d0-eta
      yp=1.d0+eta
      f(1)=.25d0*xm*ym
      f(2)=.25d0*xp*ym
      f(3)=.25d0*xp*yp
      f(4)=.25d0*xm*yp
      dfksi(1)=-.25d0*ym
      dfksi(2)= .25d0*ym
      dfksi(3)= .25d0*yp
      dfksi(4)=-.25d0*yp
      dfeta(1)=-.25d0*xm
      dfeta(2)=-.25d0*xp
      dfeta(3)= .25d0*xp
      dfeta(4)= .25d0*xm
c
         xintp(1,nn) = f(1)*x(1,1) + f(2)*x(1,2)
     *                 + f(3)*x(1,3) + f(4)*x(1,4)
         xintp(2,nn) = f(1)*x(2,1) + f(2)*x(2,2)
     *                 + f(3)*x(2,3) + f(4)*x(2,4)
         xac11 = dfksi(1)*x(1,1) + dfksi(2)*x(1,2)
     *         + dfksi(3)*x(1,3) + dfksi(4)*x(1,4)
         xac21 = dfeta(1)*x(1,1) + dfeta(2)*x(1,2)
     *         + dfeta(3)*x(1,3) + dfeta(4)*x(1,4)
         xac22 = dfeta(1)*x(2,1) + dfeta(2)*x(2,2)
     *         + dfeta(3)*x(2,3) + dfeta(4)*x(2,4)
         xac12 = dfksi(1)*x(2,1) + dfksi(2)*x(2,2)
     *         + dfksi(3)*x(2,3) + dfksi(4)*x(2,4)
         det   = xac11*xac22-xac21*xac12
       if(det.le.0.) then
         ierr=1+ierr
         return
       endif
cees det > 0.0
         dvol = det
         det = 1.d0/det
         yac11 = xac22*det
         yac22 = xac11*det
         yac21 = -xac21*det
         yac12 = -xac12*det
         dfdx(1) = yac11*dfksi(1)+yac12*dfeta(1)
         dfdx(2) = yac11*dfksi(2)+yac12*dfeta(2)
         dfdx(3) = yac11*dfksi(3)+yac12*dfeta(3)
         dfdx(4) = yac11*dfksi(4)+yac12*dfeta(4)
         dfdy(1) = yac21*dfksi(1)+yac22*dfeta(1)
         dfdy(2) = yac21*dfksi(2)+yac22*dfeta(2)
         dfdy(3) = yac21*dfksi(3)+yac22*dfeta(3)
         dfdy(4) = yac21*dfksi(4)+yac22*dfeta(4)
10    continue
c
      return
      end