Civil Engineering Reference
In-Depth Information
READ(10,*)(no(i),val(i,:),i=1,loaded_nodes)
!-----------------------load increment loop-------------------------------
READ(10,*)tol,limit,incs; ALLOCATE(qinc(incs)); READ(10,*)qinc
WRITE(11,'(/A)')
&
" step load disp iters cg iters/plastic iter"
totd=zero; tensor=zero; xnew=zero; dl=zero; diag_precon(0)=zero; ptot=zero
load_increments: DO iy=1,incs
ptot=ptot+qinc(iy); bdylds=zero; loads=zero; iters=0; cg_tot=0
DO i=1,loaded_nodes; loads(nf(:,no(i)))=val(i,:)*qinc(iy); END DO
!-----------------------plastic iteration loop --------------------------
plastic_iters: DO
iters=iters+1; IF(iters/=1)loads=zero
WRITE(*,'(A,F8.2,A,I4)')" load",ptot," iteration",iters
loads=loads+bdylds
IF(ABS(SUM(loads))<small_tol)THEN; iters=iters-1; EXIT; END IF
bdylds=zero; ddylds=zero; d=diag_precon*loads; p=d; x=zero; cg_iters=0
!-----------------------pcg equation solution-----------------------------
pcg: DO
cg_iters=cg_iters+1; u=zero
elements_3: DO iel=1,nels
g=g_g(:,iel); km=storkm(:,:,iel); u(g)=u(g)+MATMUL(km,p(g))
END DO elements_3
up=DOT_PRODUCT(loads,d); alpha=up/DOT_PRODUCT(p,u); xnew=x+p*alpha
loads=loads-u*alpha; d=diag_precon*loads
beta=DOT_PRODUCT(loads,d)/up; p=d+p*beta
call checon(xnew,x,cg_tol,cg_converged)
IF(cg_converged.OR.cg_iters==cg_limit)EXIT
END DO pcg; cg_tot=cg_tot+cg_iters; loads=xnew; loads(0)=zero
diag_precon=zero
!-----------------------go round the Gauss Points-------------------------
elements_4: DO iel=1,nels
bload=zero; dload=zero; num=g_num(:,iel)
coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); eld=loads(g); km=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)
CALL deemat(dee,prop(2,etype(iel)),prop(3,etype(iel)))
stress=tensor(:,i,iel); CALL invar(stress,sigm,dsbar,lode_theta)
ff=dsbar-SQRT(d3)*prop(1,etype(iel))
IF(ff>fftol)THEN
dlam=dl(i,iel); CALL vmflow(stress,dsbar,vmfl)
CALL fmrmat(vmfl,dsbar,dlam,dee,rmat); caflow=MATMUL(rmat,vmfl)
bot=DOT_PRODUCT(vmfl,caflow); CALL formaa(vmfl,rmat,daatd)
dee=rmat-daatd/bot
END IF; sigma=MATMUL(dee,eps); stress=sigma+tensor(:,i,iel)
CALL invar(stress,sigm,dsbar,lode_theta)
!-----------------------check whether yield is violated-------------------
fnew=dsbar-SQRT(d3)*prop(1,etype(iel)); fstiff=fnew; elso=zero
IF(fnew>=zero)THEN
CALL deemat(dee,prop(2,etype(iel)),prop(3,etype(iel)))
CALL vmflow(stress,dsbar,vmfl); caflow=MATMUL(dee,vmfl)
bot=DOT_PRODUCT(vmfl,caflow); dlam=fnew/bot; elso=caflow*dlam
stress=tensor(:,i,iel)+sigma-elso
CALL invar(stress,sigm,dsbar,lode_theta)
fnew=dsbar-SQRT(d3)*prop(1,etype(iel))
iterate_on_fnew: DO
CALL vmflow(stress,dsbar,vmfl); caflow=MATMUL(dee,vmfl)*dlam
