c***************************************************************** subroutine uelmat(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars, 1 props, nprops, coords, mcrd, nnode, u, du, v, a, jtype, time, 2 dtime, kstep, kinc, jelem, params, ndload, jdltyp, adlmag, 3 predef, npredf, lflags, mlvarx, ddlmag, mdload, pnewdt, jprops, 4 njpro, period,materiallib) c include 'aba_param.inc' C dimension rhs(mlvarx,*), amatrx(ndofel, ndofel), props(*), 1 svars(*), energy(*), coords(mcrd, nnode), u(ndofel), 2 du(mlvarx,*), v(ndofel), a(ndofel), time(2), params(*), 3 jdltyp(mdload,*), adlmag(mdload,*), ddlmag(mdload,*), 4 predef(2, npredf, nnode), lflags(*), jprops(*) parameter (zero=0.d0, one=1.d0) parameter (ndim=2, ndof=2, ndi=3, nshr=1, nnodemax=4, 1 ntens=4, ninpt=4, nsvint=8) c c ndim ... number of spatial dimensions c ndof ... number of degrees of freedom per node c ndi ... number of direct stress components c nshr ... number of shear stress component c ntens ... total number of stress tensor components c (=ndi+nshr) c ninpt ... number of integration points c nsvint... number of state variables per integration pt c (including stress and strain) c dimension stiff(ndof*nnodemax,ndof*nnodemax), 1 force(ndof*nnodemax), shape(nnodemax), dshape(ndim,nnodemax), 2 xjaci(ndim,ndim), bmat(nnodemax*ndim), statevLocal(nsvint), 3 stress(ntens), ddsdde(ntens, ntens), 4 stran(ntens), dstran(ntens), wght(ninpt) dimension coords_ip(3),dfgrd0(3,3),dfgrd1(3,3), 1 drot(3,3), ddsddt(ntens), drplde(ntens) dimension coords_new(mcrd,nnodemax) c dimension predef_loc(npredf),dpredef_loc(npredf),xx1(3,3), * xx1Old(3,3) dimension xjaci_new(ndim,ndim),bmat_new(nnodemax*ndim) c data wght /one, one, one, one/ c c******************************************************************** c c U1 = first-order, plane strain, full integration c c State variables: each integration point has nsvint SDVs c c isvinc=(npt-1)*nsvint ... integration point counter c statev(1+isvinc ) ... stress c statev(1+isvinc+ ntens) ... strain c c******************************************************************** if (lflags(3).eq.4) then do i=1, ntens do j=1, ntens ddsdde(i,j) = zero end do ddsdde(i,j) = one enddo do i=1, ndofel do j=1, ndofel amatrx(i,j) = zero end do amatrx(i,i) = one end do goto 999 end if c c Preliminaries c pnewdtLocal = pnewdt if(jtype .ne. 1) then write(7,*)'Incorrect element type' call xit endif if(nsvars .lt. ninpt*nsvint) then write(7,*)'Increase the number of SDVs to', ninpt*nsvint call xit endif c c initialize rhs and lhs c do k1=1, ndof*nnode rhs(k1, 1)= zero do k2=1, ndof*nnode amatrx(k1, k2)= zero end do end do c do k1=1,nnode do k2=1,mcrd kk = (k1-1)*mcrd + k2 coords_new(k2,k1) = coords(k2,k1) + u(kk) end do end do c c loop over integration points c do kintk = 1, ninpt c c evaluate shape functions and derivatives c call shapefcn(kintk,ninpt,nnode,ndim,shape,dshape) c c compute coordinates at the integration point c do k1=1, 3 coords_ip(k1) = zero end do do k1=1,nnode do k2=1,mcrd if (lflags(2).eq.0) then coords_ip(k2)=coords_ip(k2)+shape(k1)*coords(k2,k1) else coords_ip(k2)=coords_ip(k2)+shape(k1)*coords_new(k2,k1) end if end do end do c if(npredf.gt.0) then call tempfv(kintk,ninpt,nnode,ndim,shape,predef, * npredf,predef_loc,dpredef_loc) end if c c form B-matrix c djac = one djac_new = one call jacobian(jelem,mcrd,ndim,nnode,coords,dshape, 1 djac,xjaci,pnewdt,coords_new,xjaci_new,djac_new) if (pnewdt .lt. pnewdtLocal) pnewdtLocal = pnewdt c call bmatrix(xjaci,dshape,nnode,ndim,bmat,xjaci_new, * bmat_new) c c calculate incremental strains c call straininc(ntens,ndof,ndim,nnode,mlvarx,bmat,du,dstran, * u,xx1,xx1Old) c c call constitutive routine c call statevar(kintk,nsvint,svars,statevLocal,1) c c material routine for computing state variables c sse = zero do k1=1,ntens stress(k1) = statevLocal(k1) stran (k1) = statevLocal(k1+ntens) end do dvmat = zero celent = one call material_lib_mech(materiallib,stress,ddsdde,stran,dstran, * kintk,djac,dvmat,xx1,predef_loc,dpredef_loc,npredf,celent, * coords_ip) c do k1=1,ntens statevLocal(k1) = stress(k1) statevLocal(ntens+k1) = stran(k1) + dstran(k1) end do c call statevar(kintk,nsvint,svars,statevLocal,0) c c form stiffness matrix and internal force vector c if(lflags(2).eq.0) then call stiffmatrix(ntens,nnode,ndim,ndof, 1 wght(kintk),djac,ddsdde,stress,bmat,stiff,force) else call stiffmatrix(ntens,nnode,ndim,ndof, 1 wght(kintk),djac_new,ddsdde,stress,bmat_new,stiff,force) end if c c assemble rhs and lhs c do k1=1, ndof*nnode rhs(k1, 1) = rhs(k1, 1) - force(k1) do k2=1, ndof*nnode amatrx(k1, k2) = amatrx(k1, k2) + stiff(k1,k2) end do end do end do ! end loop on material integration points pnewdt = pnewdtLocal c 999 continue c return end c***************************************************************** c c shape functions c subroutine shapefcn(kintk,ninpt,nnode,ndim,dN,dNdz) c include 'aba_param.inc' c parameter (dmone=-1.0d0,one=1.0d0,four=4.0d0,eight=8.0d0, * gaussCoord=0.577350269d0) parameter (maxElemNode=8,maxDof=3,i2d4node=24,i3d8node=38) dimension dN(*),dNdz(ndim,*),coord24(2,4),coord38(3,8) c data coord24 /dmone, dmone, 2 one, dmone, 3 one, one, 4 dmone, one/ c data coord38 /dmone, dmone, dmone, 2 one, dmone, dmone, 3 one, one, dmone, 4 dmone, one, dmone, 5 dmone, dmone, one, 6 one, dmone, one, 7 one, one, one, 8 dmone, one, one/ C iCode = 0 if (ninpt.eq.4.and.nnode.eq.4.and.ndim.eq.2) then iCode = 24 else if (ninpt.eq.8.and.nnode.eq.8.and.ndim.eq.3) then iCode = 38 else write (6,*) '***ERROR: The shape fuctions cannot be found' end if C C 3D 8-nodes C if (iCode.eq.i3d8node) then c c determine (g,h,r) c g = coord38(1,kintk)*gaussCoord h = coord38(2,kintk)*gaussCoord r = coord38(3,kintk)*gaussCoord c c shape functions dN(1) = (one - g)*(one - h)*(one - r)/eight dN(2) = (one + g)*(one - h)*(one - r)/eight dN(3) = (one + g)*(one + h)*(one - r)/eight dN(4) = (one - g)*(one + h)*(one - r)/eight dN(5) = (one - g)*(one - h)*(one + r)/eight dN(6) = (one + g)*(one - h)*(one + r)/eight dN(7) = (one + g)*(one + h)*(one + r)/eight dN(8) = (one - g)*(one + h)*(one + r)/eight c c derivative d(Ni)/d(g) dNdz(1,1) = -(one - h)*(one - r)/eight dNdz(1,2) = (one - h)*(one - r)/eight dNdz(1,3) = (one + h)*(one - r)/eight dNdz(1,4) = -(one + h)*(one - r)/eight dNdz(1,5) = -(one - h)*(one + r)/eight dNdz(1,6) = (one - h)*(one + r)/eight dNdz(1,7) = (one + h)*(one + r)/eight dNdz(1,8) = -(one + h)*(one + r)/eight c c derivative d(Ni)/d(h) dNdz(2,1) = -(one - g)*(one - r)/eight dNdz(2,2) = -(one + g)*(one - r)/eight dNdz(2,3) = (one + g)*(one - r)/eight dNdz(2,4) = (one - g)*(one - r)/eight dNdz(2,5) = -(one - g)*(one + r)/eight dNdz(2,6) = -(one + g)*(one + r)/eight dNdz(2,7) = (one + g)*(one + r)/eight dNdz(2,8) = (one - g)*(one + r)/eight c c derivative d(Ni)/d(r) dNdz(3,1) = -(one - g)*(one - h)/eight dNdz(3,2) = -(one + g)*(one - h)/eight dNdz(3,3) = -(one + g)*(one + h)/eight dNdz(3,4) = -(one - g)*(one + h)/eight dNdz(3,5) = (one - g)*(one - h)/eight dNdz(3,6) = (one + g)*(one - h)/eight dNdz(3,7) = (one + g)*(one + h)/eight dNdz(3,8) = (one - g)*(one + h)/eight C C 2D 4-nodes C else if (iCode.eq.i2d4node) then c c determine (g,h) c g = coord24(1,kintk)*gaussCoord h = coord24(2,kintk)*gaussCoord c c shape functions dN(1) = (one - g)*(one - h)/four; dN(2) = (one + g)*(one - h)/four; dN(3) = (one + g)*(one + h)/four; dN(4) = (one - g)*(one + h)/four; c c derivative d(Ni)/d(g) dNdz(1,1) = -(one - h)/four; dNdz(1,2) = (one - h)/four; dNdz(1,3) = (one + h)/four; dNdz(1,4) = -(one + h)/four; c c derivative d(Ni)/d(h) dNdz(2,1) = -(one - g)/four; dNdz(2,2) = -(one + g)/four; dNdz(2,3) = (one + g)/four; dNdz(2,4) = (one - g)/four; end if c return end c******************************************************** subroutine tempfv(kintk,ninpt,nnode,ndim,shape,predef, * npredf,predef_loc,dpredef_loc) c include 'aba_param.inc' c dimension shape(nnode),predef(2,npredf,nnode) dimension predef_loc(npredf),dpredef_loc(npredf) parameter (zero=0.d0) c do k1=1,npredf predef_loc(k1) = zero dpredef_loc(k1) = zero do k2=1,nnode predef_loc(k1) = predef_loc(k1)+ & (predef(1,k1,k2)-predef(2,k1,k2))*shape(k2) dpredef_loc(k1) = dpredef_loc(k1)+predef(2,k1,k2)*shape(k2) end do end do c return end c***************************************************************** subroutine jacobian(jelem,mcrd,ndim,nnode, 1 coords,dshape,djac,xjaci,pnewdt,coords_new,xjaci_new, 2 djac_new) c c Notation: ndim ....... element dimension c nnode ..... number of nodes c coords ..... coordinates of nodes c dshape ..... derivs of shape fcn c djac ....... determinant of Jacobian c xjaci ...... inverse of Jacobian matrix c include 'aba_param.inc' parameter(zero=0.d0, fourth=0.25d0, maxDof=3) dimension xjac(maxDof,maxDof), xjaci(ndim,*), coords(mcrd,*) dimension dshape(ndim,*),coords_new(mcrd,*) dimension xjac_new(maxDof,maxDof), xjaci_new(ndim,*) c do i = 1, ndim do j = 1, ndim xjac(i,j) = zero xjaci(i,j) = zero xjac_new(i,j) = zero xjaci_new(i,j) = zero end do end do c do inod= 1, nnode do idim = 1, ndim do jdim = 1, ndim xjac(jdim,idim) = xjac(jdim,idim) + 1 dshape(jdim,inod)*coords(idim,inod) xjac_new(jdim,idim) = xjac_new(jdim,idim) + 1 dshape(jdim,inod)*coords_new(idim,inod) end do end do end do C C ndim == 3 C if (ndim.eq.3) then djac = xjac(1,1)*xjac(2,2)*xjac(3,3) + & xjac(2,1)*xjac(3,2)*xjac(1,3) + & xjac(3,1)*xjac(2,3)*xjac(1,2) - & xjac(3,1)*xjac(2,2)*xjac(1,3) - & xjac(2,1)*xjac(1,2)*xjac(3,3) - & xjac(1,1)*xjac(2,3)*xjac(3,2) djac_new = xjac_new(1,1)*xjac_new(2,2)*xjac_new(3,3) + & xjac_new(2,1)*xjac_new(3,2)*xjac_new(1,3) + & xjac_new(3,1)*xjac_new(2,3)*xjac_new(1,2) - & xjac_new(3,1)*xjac_new(2,2)*xjac_new(1,3) - & xjac_new(2,1)*xjac_new(1,2)*xjac_new(3,3) - & xjac_new(1,1)*xjac_new(2,3)*xjac_new(3,2) if (djac .gt. zero) then ! jacobian is positive - o.k. xjaci(1,1) = (xjac(2,2)*xjac(3,3)-xjac(2,3)*xjac(3,2))/djac xjaci(1,2) = (xjac(1,3)*xjac(3,2)-xjac(1,2)*xjac(3,3))/djac xjaci(1,3) = (xjac(1,2)*xjac(2,3)-xjac(1,3)*xjac(2,2))/djac ! xjaci(2,1) = (xjac(2,3)*xjac(3,1)-xjac(2,1)*xjac(3,3))/djac xjaci(2,2) = (xjac(1,1)*xjac(3,3)-xjac(1,3)*xjac(3,1))/djac xjaci(2,3) = (xjac(1,3)*xjac(2,1)-xjac(1,1)*xjac(2,3))/djac ! xjaci(3,1) = (xjac(2,1)*xjac(3,2)-xjac(2,2)*xjac(3,1))/djac xjaci(3,2) = (xjac(1,2)*xjac(3,1)-xjac(1,1)*xjac(3,2))/djac xjaci(3,3) = (xjac(1,1)*xjac(2,2)-xjac(1,2)*xjac(2,1))/djac else ! negative or zero jacobian write(7,*)'WARNING: element',jelem,'has neg. 1 Jacobian' pnewdt = fourth endif if (djac_new .gt. zero) then ! jacobian is positive - o.k. xjaci_new(1,1) = (xjac_new(2,2)*xjac_new(3,3)- & xjac_new(2,3)*xjac_new(3,2))/djac_new xjaci_new(1,2) = (xjac_new(1,3)*xjac_new(3,2)- & xjac_new(1,2)*xjac_new(3,3))/djac_new xjaci_new(1,3) = (xjac_new(1,2)*xjac_new(2,3)- & xjac_new(1,3)*xjac_new(2,2))/djac_new ! xjaci_new(2,1) = (xjac_new(2,3)*xjac_new(3,1)- & xjac_new(2,1)*xjac_new(3,3))/djac_new xjaci_new(2,2) = (xjac_new(1,1)*xjac_new(3,3)- & xjac_new(1,3)*xjac_new(3,1))/djac_new xjaci_new(2,3) = (xjac_new(1,3)*xjac_new(2,1)- & xjac_new(1,1)*xjac_new(2,3))/djac_new ! xjaci_new(3,1) = (xjac_new(2,1)*xjac_new(3,2)- & xjac_new(2,2)*xjac_new(3,1))/djac_new xjaci_new(3,2) = (xjac_new(1,2)*xjac_new(3,1)- & xjac_new(1,1)*xjac_new(3,2))/djac_new xjaci_new(3,3) = (xjac_new(1,1)*xjac_new(2,2)- & xjac_new(1,2)*xjac_new(2,1))/djac_new else ! negative or zero jacobian write(7,*)'WARNING: element',jelem,'has neg. 1 Jacobian' pnewdt = fourth endif C C ndim == 2 C else if (ndim.eq.2) then djac = xjac(1,1)*xjac(2,2) - xjac(1,2)*xjac(2,1) djac_new = xjac_new(1,1)*xjac_new(2,2) * - xjac_new(1,2)*xjac_new(2,1) if (djac .gt. zero) then ! jacobian is positive - o.k. xjaci(1,1) = xjac(2,2)/djac xjaci(2,2) = xjac(1,1)/djac xjaci(1,2) = -xjac(1,2)/djac xjaci(2,1) = -xjac(2,1)/djac else ! negative or zero jacobian write(7,*)'WARNING: element',jelem,'has neg. Jacobian' pnewdt = fourth endif if (djac_new .gt. zero) then ! jacobian is positive - o.k. xjaci_new(1,1) = xjac_new(2,2)/djac_new xjaci_new(2,2) = xjac_new(1,1)/djac_new xjaci_new(1,2) = -xjac_new(1,2)/djac_new xjaci_new(2,1) = -xjac_new(2,1)/djac_new else ! negative or zero jacobian write(7,*)'WARNING: element',jelem,'has neg. Jacobian' pnewdt = fourth endif end if return end c***************************************************************** subroutine bmatrix(xjaci,dshape,nnode,ndim,bmat, * xjaci_new,bmat_new) c c Notation: c bmat(i) .....dN1/dx, dN1/dy, dN2/dx, dN2/dy.. c xjaci ...... inverse Jabobian matrix c dshape ......derivative of shape functions c include 'aba_param.inc' parameter (zero=0.d0) dimension bmat(*), dshape(ndim,*) dimension xjaci(ndim,*) dimension xjaci_new(ndim,*),bmat_new(*) do i = 1, nnode*ndim bmat(i) = zero bmat_new(i) = zero end do do inod = 1, nnode do ider = 1, ndim do idim = 1, ndim irow = idim + (inod - 1)*ndim bmat(irow) = bmat(irow) + 1 xjaci(idim,ider)*dshape(ider,inod) bmat_new(irow) = bmat_new(irow) + 1 xjaci_new(idim,ider)*dshape(ider,inod) end do end do end do return end c***************************************************************** subroutine straininc(ntens,ndof,ndim,nnode, 1 mlvarx,bmat,du,dstran,u,xx1,xx1Old) c c Notation: c dstran(i) incremental strain component c note: i = 1 xx direction c = 2 yy direction c = 3 zz direction c = 4 xy direction c u() - displacement c du() - increment of displacement in the last inc. c c include 'aba_param.inc' parameter(zero=0.d0, one=1.d0) dimension dstran(*), bmat(ndim,*), 1 du(mlvarx, *), xdu(3), xx1(3,*), 2 u(mlvarx, *), utmp(3), 3 utmpOld(3),xx1Old(3,*),eps(3,3),dInvFold(3,3) do i = 1, ntens dstran(i) = zero end do ! ! set xx1 to Identity matrix do k1=1,3 do k2=1,3 xx1(k1,k2) = zero xx1Old(k1,k2) = zero end do xx1(k1,k1) = one xx1Old(k1,k1) = one end do c c************************************ c Compute incremental strains c************************************ c do nodi = 1, nnode incr_row = (nodi - 1)*ndof do i = 1, ndof xdu(i)= du(i + incr_row,1) utmp(i) = u(i + incr_row,1) utmpOld(i) = utmp(i)-xdu(i) end do dNidx = bmat(1,nodi) dNidy = bmat(2,nodi) if (ndim.eq.3) then dNidz = bmat(3,nodi) dstran(1) = dstran(1) + dNidx*xdu(1) dstran(2) = dstran(2) + dNidy*xdu(2) dstran(3) = dstran(3) + dNidz*xdu(3) ! (1-2) dstran(4) = dstran(4) + 1 dNidy*xdu(1) + 2 dNidx*xdu(2) ! (1-3) dstran(5) = dstran(5) + 1 dNidz*xdu(1) + 2 dNidx*xdu(3) ! (2-3) dstran(6) = dstran(6) + 1 dNidy*xdu(3) + 2 dNidz*xdu(2) c deformation gradient xx1(1,1) = xx1(1,1) + dNidx*utmp(1) xx1(1,2) = xx1(1,2) + dNidy*utmp(1) xx1(1,3) = xx1(1,3) + dNidz*utmp(1) ! xx1(2,1) = xx1(2,1) + dNidx*utmp(2) xx1(2,2) = xx1(2,2) + dNidy*utmp(2) xx1(2,3) = xx1(2,3) + dNidz*utmp(2) ! xx1(3,1) = xx1(3,1) + dNidx*utmp(3) xx1(3,2) = xx1(3,2) + dNidy*utmp(3) xx1(3,3) = xx1(3,3) + dNidz*utmp(3) else if (ndim.eq.2) then dstran(1) = dstran(1) + dNidx*xdu(1) dstran(2) = dstran(2) + dNidy*xdu(2) dstran(4) = dstran(4) + 1 dNidy*xdu(1) + 2 dNidx*xdu(2) c deformation gradient xx1(1,1) = xx1(1,1) + dNidx*utmp(1) xx1(1,2) = xx1(1,2) + dNidy*utmp(1) xx1(2,1) = xx1(2,1) + dNidx*utmp(2) xx1(2,2) = xx1(2,2) + dNidy*utmp(2) c xx1Old(1,1) = xx1Old(1,1) + dNidx*utmpOld(1) xx1Old(1,2) = xx1Old(1,2) + dNidy*utmpOld(1) xx1Old(2,1) = xx1Old(2,1) + dNidx*utmpOld(2) xx1Old(2,2) = xx1Old(2,2) + dNidy*utmpOld(2) end if end do c return end c***************************************************************** subroutine statevar(npt,nsvint,statev,statev_ip,icopy) c c Transfer data to/from element-level state variable array from/to c material-point level state variable array. c include 'aba_param.inc' dimension statev(*),statev_ip(*) isvinc= (npt-1)*nsvint ! integration point increment if (icopy .eq. 1) then c c Prepare arrays for entry into umat c do i = 1, nsvint statev_ip(i)=statev(i+isvinc) end do c else c c Update element state variables upon return from umat c do i = 1, nsvint statev(i+isvinc)=statev_ip(i) end do end if return end c***************************************************************** subroutine stiffmatrix(ntens,nnode,ndim,ndof, 1 weight,djac,ddsdde,stress,bmat,stiff,force) c c Stiffness matrix and internal force contributions at c material integration point c include 'aba_param.inc' parameter(zero=0.d0,maxDof=3) dimension stiff(ndof*nnode,*), stiff_p(maxDof,maxDof) dimension force(*), force_p(maxDof) dimension stress(*),bmat(ndim,*),ddsdde(ntens,*) dNjdx = zero dNjdy = zero dNjdz = zero do i = 1, ndof*nnode force(i) = zero do j = 1, ndof*nnode stiff(j,i) = zero end do end do dvol= weight*djac do nodj = 1, nnode incr_col = (nodj - 1)*ndof dNjdx = bmat(1,nodj) dNjdy = bmat(2,nodj) if (ndim.eq.2) then force_p(1) = dNjdx*stress(1) + dNjdy*stress(4) force_p(2) = dNjdy*stress(2) + dNjdx*stress(4) else if (ndim.eq.3) then dNjdz = bmat(3,nodj) force_p(1) = dNjdx*stress(1) + & dNjdy*stress(4) + & dNjdz*stress(5) force_p(2) = dNjdy*stress(2) + & dNjdx*stress(4) + & dNjdz*stress(6) force_p(3) = dNjdy*stress(6) + & dNjdx*stress(5) + & dNjdz*stress(3) end if do jdof = 1, ndof jcol = jdof + incr_col force(jcol) = force(jcol) + & force_p(jdof)*dvol end do do nodi = 1, nnode incr_row = (nodi -1)*ndof dNidx = bmat(1,nodi) dNidy = bmat(2,nodi) if (ndim.eq.2) then stiff_p(1,1) = dNidx*ddsdde(1,1)*dNjdx & + dNidy*ddsdde(4,4)*dNjdy & + dNidx*ddsdde(1,4)*dNjdy & + dNidy*ddsdde(4,1)*dNjdx stiff_p(1,2) = dNidx*ddsdde(1,2)*dNjdy & + dNidy*ddsdde(4,4)*dNjdx & + dNidx*ddsdde(1,4)*dNjdx & + dNidy*ddsdde(4,2)*dNjdy stiff_p(2,1) = dNidy*ddsdde(2,1)*dNjdx & + dNidx*ddsdde(4,4)*dNjdy & + dNidy*ddsdde(2,4)*dNjdy & + dNidx*ddsdde(4,1)*dNjdx stiff_p(2,2) = dNidy*ddsdde(2,2)*dNjdy & + dNidx*ddsdde(4,4)*dNjdx & + dNidy*ddsdde(2,4)*dNjdx & + dNidx*ddsdde(4,2)*dNjdy else if (ndim.eq.3) then dNidz = bmat(3,nodi) stiff_p(1,1) = dNjdx*ddsdde(1,1)*dNidx & + dNjdx*ddsdde(1,4)*dNidy & + dNjdx*ddsdde(1,5)*dNidz & + dNjdy*ddsdde(4,1)*dNidx & + dNjdy*ddsdde(4,4)*dNidy & + dNjdy*ddsdde(4,5)*dNidz & + dNjdz*ddsdde(5,1)*dNidx & + dNjdz*ddsdde(5,4)*dNidy & + dNjdz*ddsdde(5,5)*dNidz c stiff_p(1,2) = dNjdx*ddsdde(1,2)*dNidy & + dNjdx*ddsdde(1,4)*dNidx & + dNjdx*ddsdde(1,6)*dNidz & + dNjdy*ddsdde(4,2)*dNidy & + dNjdy*ddsdde(4,4)*dNidx & + dNjdy*ddsdde(4,6)*dNidz & + dNjdz*ddsdde(5,2)*dNidy & + dNjdz*ddsdde(5,4)*dNidx & + dNjdz*ddsdde(5,6)*dNidz c stiff_p(1,3) = dNjdx*ddsdde(1,3)*dNidz & + dNjdx*ddsdde(1,5)*dNidx & + dNjdx*ddsdde(1,6)*dNidy & + dNjdy*ddsdde(4,3)*dNidz & + dNjdy*ddsdde(4,5)*dNidx & + dNjdy*ddsdde(4,6)*dNidy & + dNjdz*ddsdde(5,3)*dNidz & + dNjdz*ddsdde(5,5)*dNidx & + dNjdz*ddsdde(5,6)*dNidy c stiff_p(2,1) = dNjdy*ddsdde(2,1)*dNidx & + dNjdy*ddsdde(2,4)*dNidy & + dNjdy*ddsdde(2,5)*dNidz & + dNjdx*ddsdde(4,1)*dNidx & + dNjdx*ddsdde(4,4)*dNidy & + dNjdx*ddsdde(4,5)*dNidz & + dNjdz*ddsdde(6,1)*dNidx & + dNjdz*ddsdde(6,4)*dNidy & + dNjdz*ddsdde(6,5)*dNidz c stiff_p(2,2) = dNjdy*ddsdde(2,2)*dNidy & + dNjdy*ddsdde(2,4)*dNidx & + dNjdy*ddsdde(2,6)*dNidz & + dNjdx*ddsdde(4,2)*dNidy & + dNjdx*ddsdde(4,4)*dNidx & + dNjdx*ddsdde(4,6)*dNidz & + dNjdz*ddsdde(6,2)*dNidy & + dNjdz*ddsdde(6,4)*dNidx & + dNjdz*ddsdde(6,6)*dNidz c stiff_p(2,3) = dNjdy*ddsdde(2,3)*dNidz & + dNjdy*ddsdde(2,5)*dNidx & + dNjdy*ddsdde(2,6)*dNidy & + dNjdx*ddsdde(4,3)*dNidz & + dNjdx*ddsdde(4,5)*dNidx & + dNjdx*ddsdde(4,6)*dNidy & + dNjdz*ddsdde(6,3)*dNidz & + dNjdz*ddsdde(6,5)*dNidx & + dNjdz*ddsdde(6,6)*dNidy c stiff_p(3,1) = dNjdz*ddsdde(3,1)*dNidx & + dNjdz*ddsdde(3,4)*dNidy & + dNjdz*ddsdde(3,5)*dNidz & + dNjdx*ddsdde(5,1)*dNidx & + dNjdx*ddsdde(5,4)*dNidy & + dNjdx*ddsdde(5,5)*dNidz & + dNjdy*ddsdde(6,1)*dNidx & + dNjdy*ddsdde(6,4)*dNidy & + dNjdy*ddsdde(6,5)*dNidz c stiff_p(3,2) = dNjdz*ddsdde(3,2)*dNidy & + dNjdz*ddsdde(3,4)*dNidx & + dNjdz*ddsdde(3,6)*dNidz & + dNjdx*ddsdde(5,2)*dNidy & + dNjdx*ddsdde(5,4)*dNidx & + dNjdx*ddsdde(5,6)*dNidz & + dNjdy*ddsdde(6,2)*dNidy & + dNjdy*ddsdde(6,4)*dNidx & + dNjdy*ddsdde(6,6)*dNidz c stiff_p(3,3) = dNjdz*ddsdde(3,3)*dNidz & + dNjdz*ddsdde(3,5)*dNidx & + dNjdz*ddsdde(3,6)*dNidy & + dNjdx*ddsdde(5,3)*dNidz & + dNjdx*ddsdde(5,5)*dNidx & + dNjdx*ddsdde(5,6)*dNidy & + dNjdy*ddsdde(6,3)*dNidz & + dNjdy*ddsdde(6,5)*dNidx & + dNjdy*ddsdde(6,6)*dNidy end if do jdof = 1, ndof icol = jdof + incr_col do idof = 1, ndof irow = idof + incr_row stiff(irow,icol) = stiff(irow,icol) + & stiff_p(idof,jdof)*dvol end do end do end do end do c return end