Civil Engineering Reference
In-Depth Information
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
CALL emb_3d_bc(ifix,nx1,nx2,ny1,ny2,nze,nf); neq=MAXVAL(nf)
ALLOCATE(loads(0:neq),bdylds(0:neq),oldis(0:neq),gravlo(0:neq),kdiag(neq))
!-----------------------loop the elements to find global arrays sizes-----
kdiag=0
elements_1: DO iel=1,nels
CALL emb_3d_geom(iel,nx1,nx2,ny1,ny2,nze,w1,s1,w2,h1,h2,d1,coord,num)
CALL num_to_g(num,nf,g); g_num(:,iel)=num
g_coord(:,num)=TRANSPOSE(coord); g_g(:,iel)=g; CALL fkdiag(kdiag,g)
END DO elements_1
DO i=2,neq; kdiag(i)=kdiag(i)+kdiag(i-1); END DO; ALLOCATE(kv(kdiag(neq)))
WRITE(11,'(A,I7,A,I8)')
&
" There are",neq," equations and the skyline storage is",kdiag(neq)
CALL sample(element,points,weights); kv=zero; gravlo=zero
!-----------------------element stiffness integration and assembly--------
elements_2: DO iel=1,nels
CALL deemat(dee,prop(5,etype(iel)),prop(6,etype(iel))); num=g_num(:,iel)
coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); km=zero; eld=zero
gauss_pts_1: DO i=1,nip
CALL shape_fun(fun,points,i); 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)
eld(nodof:ndof:nodof)=eld(nodof:ndof:nodof)+fun(:)*det*weights(i)
END DO gauss_pts_1
CALL fsparv(kv,km,g,kdiag); gravlo(g)=gravlo(g)-eld*prop(4,etype(iel))
END DO elements_2
!-----------------------factorise equations-------------------------------
CALL sparin(kv,kdiag); pi=ACOS(-one)
!-----------------------trial strength reduction factor loop--------------
READ(10,*)tol,limit,nsrf; ALLOCATE(srf(nsrf)); READ(10,*)srf
WRITE(11,'(/A)')"
srf
max disp iters"
srf_trials: DO iy=1,nsrf
dt=start_dt
DO i=1,np_types
phi=prop(1,i); tnph=TAN(phi*pi/d180); phif=ATAN(tnph/srf(iy))
snph=SIN(phif); e=prop(5,i); v=prop(6,i)
ddt=d4*(one+v)*(one-two*v)/(e*(one-two*v+snph**2)); IF(ddt<dt)dt=ddt
END DO; iters=0; bdylds=zero; evpt=zero; oldis=zero
!-----------------------plastic iteration loop----------------------------
its: DO
fmax=zero; iters=iters+1; loads=gravlo+bdylds
CALL spabac(kv,loads,kdiag); loads(0)=zero
!-----------------------check plastic 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_3: DO iel=1,nels
bload=zero; phi=prop(1,etype(iel)); tnph=TAN(phi*pi/d180)
phif=ATAN(tnph/srf(iy))*d180/pi; psi=prop(3,etype(iel))
tnps=TAN(psi*pi/d180); psif=ATAN(tnps/srf(iy))*d180/pi
cf=prop(2,etype(iel))/srf(iy); e=prop(5,etype(iel))
v=prop(6,etype(iel)); CALL deemat(dee,e,v); num=g_num(:,iel)
coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); eld=loads(g)
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)
