Coupled Problems - Programming the Finite Element Method - page 426

!-----------------------check plastic convergence-------------------------
newdis=loads; newdis(nf(3,:))=zero
CALL checon(newdis,oldis,tol,converged); IF(iters==1)converged=.FALSE.
IF(converged.OR.iters==limit)bdylds=zero
!-----------------------go round the Gauss Points ------------------------
elements_4: DO iel=1,nels
phi=prop(5,etype(iel)); coh=prop(6,etype(iel))
psi=prop(7,etype(iel))
CALL deemat(dee,prop(3,etype(iel)),prop(4,etype(iel)))
num=g_num(:,iel); coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel)
eld=loads(g(1:ndof)); bload=zero
gauss_points_2: DO i=1,nip
CALL shape_der(der,points,i)
jac=MATMUL(der,coord); det=determinant(jac); CALL invert(jac)
deriv=MATMUL(jac,der); CALL beemat(bee,deriv)
eps=MATMUL(bee,eld); eps=eps-evpt(:,i,iel); sigma=MATMUL(dee,eps)
stress=sigma+tensor(1:4,i,iel)
CALL invar(stress,sigm,dsbar,lode_theta)
!-----------------------check whether yield is violated-------------------
CALL mocouf(phi,coh,sigm,dsbar,lode_theta,f)
IF(converged.OR.iters==limit)THEN; devp=stress; ELSE
IF(f>=zero)THEN
CALL mocouq(psi,dsbar,lode_theta,dq1,dq2,dq3)
CALL formm(stress,m1,m2,m3)
flow=f*(m1*dq1+m2*dq2+m3*dq3); erate=MATMUL(flow,stress)
evp=erate*dt; evpt(:,i,iel)=evpt(:,i,iel)+evp
devp=MATMUL(dee,evp)
END IF
END IF
IF(f>=zero)THEN; eload=MATMUL(TRANSPOSE(bee),devp)
bload=bload+eload*det*weights(i)
END IF
IF(converged.OR.iters==limit)THEN
!--------------update the Gauss Point stresses and porepressures----------
tensor(1:4,i,iel)=stress; dpore=zero
CALL shape_fun(funf,points,i)
DO k=1,nodf; dpore=dpore+funf(k)*loads(g(k+ndof)); END DO
tensor(5,i,iel)=tensor(5,i,iel)+dpore
END IF
END DO gauss_points_2
!-----------------------compute the total bodyloads vector ---------------
bdylds(g(1:ndof))=bdylds(g(1:ndof))+bload; bdylds(0)=zero
END DO elements_4; IF(converged.OR.iters==limit)EXIT
END DO its; disps=disps+loads
IF(j/npri*npri==j.OR.iters==limit)WRITE(11,'(5E12.4,I5)') &

time,tot_load,disps(nf(:,nres)),iters
IF(iters==limit)EXIT
END DO time_steps
CALL dismsh(loads,nf,0.05_iwp,g_coord,g_num,13)
CALL vecmsh(loads,nf,0.05_iwp,0.1_iwp,g_coord,g_num,14)
STOP
END PROGRAM p94

New scalar reals:
coh soil cohesion
cons consolidating stress ( σ 3 )
ddt

used to find the critical time step