Civil Engineering Reference
In-Depth Information
tensor(1,i,iel)=k0*tensor(2,i,iel)
tensor(4,i,iel)=tensor(1,i,iel); tensor(3,i,iel)=zero
END IF
CALL shape_der(der,points,i)
jac=MATMUL(der,coord); det=determinant(jac); CALL invert(jac)
deriv=MATMUL(jac,der); CALL beemat(bee,deriv)
km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
DO k=2,ndof,2; eld(k)=eld(k)+fun(k/2)*det*weights(i); END DO
END DO gauss_pts_1
CALL fsparv(kv,km,g,kdiag)
IF(ii<=lifts)gravlo(g)=gravlo(g)-eld*gamma; gravlo(0)=zero
END DO elements_3
!-----------------------factorise equations--and-factor gravlo by incs----
CALL sparin(kv,kdiag); gravlo=gravlo/incs
!-----------------------apply gravity loads incrementally-----------------
load_incs: DO iy=1,incs
iters=0; oldis=zero; bdylds=zero; evpt(:,:,1:oldele)=zero
!-----------------------iteration loop------------------------------------
its: DO
iters=iters+1
WRITE(*,'(A,I3,A,I3,A,I4)')" lift",ii," increment",iy,
&
" iteration",iters
loads=zero; loads=gravlo+bdylds; CALL spabac(kv,loads,kdiag)
!-----------------------check convergence---------------------------------
CALL checon(loads,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,oldele
IF(etype(iel)==1)THEN; phi=phi_f; c=c_f; e=e_f; v=v_f; psi=psi_f
ELSE; phi=phi_e; c=c_e; e=e_e; v=v_e; psi=psi_e
END IF; snph=SIN(phi*pi/d180)
dt=d4*(one+v)*(one-two*v)/(e*(one-two*v+snph**2))
CALL deemat(dee,e,v); bload=zero; num=g_num(:,iel)
coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); eld=loads(g)
gauss_pts_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)
IF(ii==1)THEN; stress=tensor(:,i,iel)
ELSE; stress=tensor(:,i,iel)+sigma
END IF; CALL invar(stress,sigm,dsbar,lode_theta)
!-----------------------check whether yield is violated-------------------
CALL mocouf(phi,c,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(devp,bee)
bload=bload+eload*det*weights(i)
END IF
!-----------------------if appropriate update the Gauss point stresses----
