Hi everyone,

I'm trying to delete elements from my mesh when a critical strain is reached. To be able to do that I wrote the following subroutine:

subroutine vusdfld(nblock,nstatev,nfieldv,nprops,ndir, & nshr,jElem,kIntPt,kLayer,kSecPt, & stepTime,totalTime,dt,cmname, & coordMp,direct,T,charLength,props, & stateOld,stateNew,field ) C C include 'ABA_PARAM.INC' C include 'vaba_param.inc' C dimension jElem(nblock), coordMp(nblock,*), 1 direct(nblock,3,3), T(nblock,3,3), 2 charLength(nblock), props(nprops), 3 stateOld(nblock,nstatev), 4 stateNew(nblock,nstatev), 5 field(nblock,nfieldv) character*80 cmname

parameter( nrData=6 ) character*3 cData(maxblk*nrData) dimension rData(maxblk*nrData), jData(maxblk*nrData) C real*8,dimension(6)::sigmaNew real*8,dimension(1)::Critical_Strain real*8,dimension(1)::MaxStress real*8,dimension(1)::i1 real*8,dimension(1)::i2 real*8,dimension(1)::i3 real*8,dimension(1)::term1 real*8,dimension(1)::term2 real*8,dimension(1)::term3 real*8,dimension(1)::termaux real*8,dimension(1)::phiaux real*8,dimension(1)::phi real*8,dimension(1)::maxprincipalE C jStatus = 1 call vgetvrm('LE',rdata,jdata,cdata,jStatus) if( jStatus .ne. 0 ) then call xplb_abqerr(-2,'Utility routine VGETVRM '// . 'failed to get variable.',0,zero,' ') call xplb_exit end if C Critical_Strain(1)=0.3000 C do k = 1, nblock StrainE11 = rData(k) StrainE22 = rData(nblock+k) StrainE33 = rData(2*nblock+k) StrainE12 = rData(3*nblock+k) StrainE23 = rData(4*nblock+k) StrainE31 = rData(5*nblock+k)

i1(1) = StrainE11 + StrainE22 + StrainE33 i2(1) = StrainE11*StrainE22 + StrainE22*StrainE33 + StrainE33*StrainE11 - StrainE12*StrainE12 - StrainE23*StrainE23 - StrainE31*StrainE31 i3(1) = StrainE11*StrainE22*StrainE33 - StrainE11*StrainE23*StrainE23 - StrainE22*StrainE31*StrainE31 - StrainE33*StrainE12*StrainE12 + 2*StrainE12*StrainE23*StrainE31

phiaux(1) = (2*i1(1)*i1(1)*i1(1) - 9*i1(1)*i2(1) + 27*i3(1)*i3(1)*i3(1))/(2*(i1(1)*i1(1)-3*i2(1))**(3./2.)) phi(1) = 1/3 * 1/cos(phiaux(1))

termaux(1) = i1(1)*i1(1) - 3*i2(1)*i2(1) term1(1) = i1(1)/3 + 2/3*SQRT(termaux(1))*cos(phi(1)) term2(1) = i1(1)/3 + 2/3*SQRT(termaux(1))*cos(phi(1) - 2*3.141592/3.) term3(1) = i1(1)/3 + 2/3*SQRT(termaux(1))*cos(phi(1) - 4*3.141592/3.)

maxprincipalE(1) = term1(1) + term2(1) + term3(1)

if (maxprincipalE(1) .ge. Critical_Strain(1)) then stateNew(k,1) = 0 end if end do RETURN END

Do you see any mistake in that? I can't have it compiled and also can't find the error.

Thanks very much!