C tot_energy_3d_1.f subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,ni,layer,nndi, * nnshear,jpltcd) C C This subroutine computes the volume and total strain energy of a C 8 noded 3D element mesh. C Post code -1 Volume C Post code -2 Strain energy density C Post code -3 Total strain energy c* * * * * * c c select a variable contour plotting (user subroutine). c c v variable c s (idss) stress array c sp stresses in preferred direction c etot total strain (generalized) c eplas total plastic strain c ecreep total creep strain c t current temperature c m element number c nn integration point number c layer layer number c ndi (3) number of direct stress components c nshear (3) number of shear stress components c c* * * * * * implicit real*8 (a-h,o-z) dp dimension s(1),etot(1),eplas(1),ecreep(1),sp(1) dimension xintp(3,1),f(8),dfdx(8),dfdy(8),dfdz(8),dvol(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 increment",incold," is",total total=0. write(75,*) endif incold = inc C Obtain strain energy density C ********************************** lms = nnx-1 lms2 = lms*neqst mms2 = lms2 lofr = (m-1)*nelstr jtoten = itoten + lofr + mms2 totde = vars(jtoten) C ********************************** C Compute element volume la1=icoord call ppp007(f,dfdx,dfdy,dfdz,vars(la1),dvol,xintp,m,ni,ierr) vol=8.*dvol(1) if(jpltcd .eq. 1)then v=vol elseif (jpltcd .eq. 2)then v=totde elseif (jpltcd .eq. 3)then v=vol*totde total = total+v write(76,*)"Increment",inc," Total Energy",total endif return end C subroutine ppp007(f,dfdx,dfdy,dfdz,x,dvol,xintp,m,nn,ierr) implicit real*8 (a-h,o-z) c dimension x(3,1),xintp(3,1),f(8),dfdx(8),dfdy(8),dfdz(8), $ gausx(8),gausy(8),gausz(8),dfksi(8),dfeta(8), $ dfzeta(8),dvol(1) c data gausx/-.577350269, .577350269,-.577350269, .577350269, $-.577350269, .577350269,-.577350269, .577350269/ data gausy/-.577350269,-.577350269, .577350269, .577350269, $-.577350269,-.577350269, .577350269, .577350269/ data gausz/-.577350269,-.577350269, -.577350269, -.577350269, $+.577350269,+.577350269, .577350269, .577350269/ c xsi = gausx(nn) eta = gausy(nn) zeta = gausz(nn) xm=1.0d0-xsi xp=1.0d0+xsi ym=1.0d0-eta yp=1.0d0+eta zm=1.d0-zeta zp=1.d0+zeta f(1)=0.125d0*xm*ym*zm f(2)=0.125d0*xp*ym*zm f(3)=0.125d0*xp*yp*zm f(4)=0.125d0*xm*yp*zm f(5)=0.125d0*xm*ym*zp f(6)=0.125d0*xp*ym*zp f(7)=0.125d0*xp*yp*zp f(8)=0.125d0*xm*yp*zp dfksi(1)=-0.125d0*ym*zm dfksi(2)= 0.125d0*ym*zm dfksi(3)= 0.125d0*yp*zm dfksi(4)=-0.125d0*yp*zm dfksi(5)=-0.125d0*ym*zp dfksi(6)= 0.125d0*ym*zp dfksi(7)= 0.125d0*yp*zp dfksi(8)=-0.125d0*yp*zp dfeta(1)=-0.125d0*xm*zm dfeta(2)=-0.125d0*xp*zm dfeta(3)= 0.125d0*xp*zm dfeta(4)= 0.125d0*xm*zm dfeta(5)=-0.125d0*xm*zp dfeta(6)=-0.125d0*xp*zp dfeta(7)= 0.125d0*xp*zp dfeta(8)= 0.125d0*xm*zp dfzeta(1) = -0.125d0*xm*ym dfzeta(2) = -0.125d0*xp*ym dfzeta(3) = -0.125d0*xp*yp dfzeta(4) = -0.125d0*xm*yp dfzeta(5) = 0.125d0*xm*ym dfzeta(6) = 0.125d0*xp*ym dfzeta(7) = 0.125d0*xp*yp dfzeta(8) = 0.125d0*xm*yp c xintp(1,nn) = f(1)*x(1,1) + f(2)*x(1,2) * + f(3)*x(1,3) + f(4)*x(1,4) * + f(5)*x(1,5) + f(6)*x(1,6) * + f(7)*x(1,7) + f(8)*x(1,8) xintp(2,nn) = f(1)*x(2,1) + f(2)*x(2,2) * + f(3)*x(2,3) + f(4)*x(2,4) * + f(5)*x(2,5) + f(6)*x(2,6) * + f(7)*x(2,7) + f(8)*x(2,8) xintp(3,nn) = f(1)*x(3,1) + f(2)*x(3,2) * + f(3)*x(3,3) + f(4)*x(3,4) * + f(5)*x(3,5) + f(6)*x(3,6) * + f(7)*x(3,7) + f(8)*x(3,8) xac11 = dfksi(1)*x(1,1) + dfksi(2)*x(1,2) * + dfksi(3)*x(1,3) + dfksi(4)*x(1,4) * + dfksi(5)*x(1,5) + dfksi(6)*x(1,6) * + dfksi(7)*x(1,7) + dfksi(8)*x(1,8) xac21 = dfeta(1)*x(1,1) + dfeta(2)*x(1,2) * + dfeta(3)*x(1,3) + dfeta(4)*x(1,4) * + dfeta(5)*x(1,5) + dfeta(6)*x(1,6) * + dfeta(7)*x(1,7) + dfeta(8)*x(1,8) xac31 = dfzeta(1)*x(1,1) + dfzeta(2)*x(1,2) * + dfzeta(3)*x(1,3) + dfzeta(4)*x(1,4) * + dfzeta(5)*x(1,5) + dfzeta(6)*x(1,6) * + dfzeta(7)*x(1,7) + dfzeta(8)*x(1,8) xac22 = dfeta(1)*x(2,1) + dfeta(2)*x(2,2) * + dfeta(3)*x(2,3) + dfeta(4)*x(2,4) * + dfeta(5)*x(2,5) + dfeta(6)*x(2,6) * + dfeta(7)*x(2,7) + dfeta(8)*x(2,8) xac12 = dfksi(1)*x(2,1) + dfksi(2)*x(2,2) * + dfksi(3)*x(2,3) + dfksi(4)*x(2,4) * + dfksi(5)*x(2,5) + dfksi(6)*x(2,6) * + dfksi(7)*x(2,7) + dfksi(8)*x(2,8) xac32 = dfzeta(1)*x(2,1) + dfzeta(2)*x(2,2) * + dfzeta(3)*x(2,3) + dfzeta(4)*x(2,4) * + dfzeta(5)*x(2,5) + dfzeta(6)*x(2,6) * + dfzeta(7)*x(2,7) + dfzeta(8)*x(2,8) xac23 = dfeta(1)*x(3,1) + dfeta(2)*x(3,2) * + dfeta(3)*x(3,3) + dfeta(4)*x(3,4) * + dfeta(5)*x(3,5) + dfeta(6)*x(3,6) * + dfeta(7)*x(3,7) + dfeta(8)*x(3,8) xac13 = dfksi(1)*x(3,1) + dfksi(2)*x(3,2) * + dfksi(3)*x(3,3) + dfksi(4)*x(3,4) * + dfksi(5)*x(3,5) + dfksi(6)*x(3,6) * + dfksi(7)*x(3,7) + dfksi(8)*x(3,8) xac33 = dfzeta(1)*x(3,1) + dfzeta(2)*x(3,2) * + dfzeta(3)*x(3,3) + dfzeta(4)*x(3,4) * + dfzeta(5)*x(3,5) + dfzeta(6)*x(3,6) * + dfzeta(7)*x(3,7) + dfzeta(8)*x(3,8) c cof11 = xac22*xac33-xac32*xac23 cof12 = xac32*xac13-xac12*xac33 cof13 = xac12*xac23-xac22*xac13 cof21 = xac23*xac31-xac33*xac21 cof22 = xac33*xac11-xac13*xac31 cof23 = xac13*xac21-xac23*xac11 cof31 = xac21*xac32-xac31*xac22 cof32 = xac31*xac12-xac11*xac32 cof33 = xac11*xac22-xac21*xac12 det = xac11*cof11 + xac21*cof12 + xac31*cof13 cees det > 0.0 if(det.le.0.) then ierr=1+ierr return endif dvol(nn) = det det = 1.d0/det dfdx(1) = (cof11*dfksi(1)+cof12*dfeta(1)+cof13*dfzeta(1))*det dfdx(2) = (cof11*dfksi(2)+cof12*dfeta(2)+cof13*dfzeta(2))*det dfdx(3) = (cof11*dfksi(3)+cof12*dfeta(3)+cof13*dfzeta(3))*det dfdx(4) = (cof11*dfksi(4)+cof12*dfeta(4)+cof13*dfzeta(4))*det dfdx(5) = (cof11*dfksi(5)+cof12*dfeta(5)+cof13*dfzeta(5))*det dfdx(6) = (cof11*dfksi(6)+cof12*dfeta(6)+cof13*dfzeta(6))*det dfdx(7) = (cof11*dfksi(7)+cof12*dfeta(7)+cof13*dfzeta(7))*det dfdx(8) = (cof11*dfksi(8)+cof12*dfeta(8)+cof13*dfzeta(8))*det dfdy(1) = (cof21*dfksi(1)+cof22*dfeta(1)+cof23*dfzeta(1))*det dfdy(2) = (cof21*dfksi(2)+cof22*dfeta(2)+cof23*dfzeta(2))*det dfdy(3) = (cof21*dfksi(3)+cof22*dfeta(3)+cof23*dfzeta(3))*det dfdy(4) = (cof21*dfksi(4)+cof22*dfeta(4)+cof23*dfzeta(4))*det dfdy(5) = (cof21*dfksi(5)+cof22*dfeta(5)+cof23*dfzeta(5))*det dfdy(6) = (cof21*dfksi(6)+cof22*dfeta(6)+cof23*dfzeta(6))*det dfdy(7) = (cof21*dfksi(7)+cof22*dfeta(7)+cof23*dfzeta(7))*det dfdy(8) = (cof21*dfksi(8)+cof22*dfeta(8)+cof23*dfzeta(8))*det dfdz(1) = (cof31*dfksi(1)+cof32*dfeta(1)+cof33*dfzeta(1))*det dfdz(2) = (cof31*dfksi(2)+cof32*dfeta(2)+cof33*dfzeta(2))*det dfdz(3) = (cof31*dfksi(3)+cof32*dfeta(3)+cof33*dfzeta(3))*det dfdz(4) = (cof31*dfksi(4)+cof32*dfeta(4)+cof33*dfzeta(4))*det dfdz(5) = (cof31*dfksi(5)+cof32*dfeta(5)+cof33*dfzeta(5))*det dfdz(6) = (cof31*dfksi(6)+cof32*dfeta(6)+cof33*dfzeta(6))*det dfdz(7) = (cof31*dfksi(7)+cof32*dfeta(7)+cof33*dfzeta(7))*det dfdz(8) = (cof31*dfksi(8)+cof32*dfeta(8)+cof33*dfzeta(8))*det 10 continue c return end