Civil Engineering Reference
In-Depth Information
kc=kc+MATMUL(MATMUL(TRANSPOSE(derivf),kay),derivf)*det*weights(i)*dtim
CALL cross_product(vol,funf,volf); c=c+volf*det*weights(i)
END DO gauss_points_1; CALL fmkdke(km,kc,c,ke,kd,theta)
storke(:,:,iel)=ke; storkd(:,:,iel)=kd
DO k=1,ndof; diag_precon(g(k))=diag_precon(g(k))+theta*km(k,k); END DO
DO k=1,nodf
diag_precon(g(ndof+k))=diag_precon(g(ndof+k))-theta*theta*kc(k,k)
END DO
END DO elements_2
diag_precon(1:)=one/diag_precon(1:); diag_precon(0)=zero
!-----------------------time stepping loop--------------------------------
WRITE(11,'(/A,I5)')" Results at node",nres; loads=zero
WRITE(11,'(4X,A)')
&
"time
load
x-disp
y-disp
porepressure cg iters"
WRITE(11,'(5E12.4)')0.0,0.0,loads(nf(:,nres))
tot_load=zero; time=zero
time_steps: DO j=1,nstep
ans=zero; time=time+dtim
elements_3: DO iel=1,nels
g=g_g(:,iel); kd=storkd(:,:,iel); ans(g)=ans(g)+MATMUL(kd,loads(g))
END DO elements_3; ans(0)=zero
!-----------------------apply absolute loading----------------------------
DO i=1,loaded_nodes
ans(nf(1:2,no(i)))=ans(nf(1:2,no(i)))+val(i,:)*(tot_load+theta*al(j))
END DO
tot_load=tot_load+al(j); d=diag_precon*ans; p=d; x=zero; cg_iters=0
!-----------------------pcg equation solution-----------------------------
pcg: DO
cg_iters=cg_iters+1; u=zero
elements_4: DO iel=1,nels
g=g_g(:,iel); ke=storke(:,:,iel); u(g)=u(g)+MATMUL(ke,p(g))
END DO elements_4
up=DOT_PRODUCT(ans,d); alpha=up/DOT_PRODUCT(p,u); xnew=x+p*alpha
ans=ans-u*alpha; d=diag_precon*ans; beta=DOT_PRODUCT(ans,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; loads=xnew; loads(0)=zero
IF(j/npri*npri==j)WRITE(11,'(5E12.4,I7)')
&
time,tot_load,loads(nf(:,nres)),cg_iters
!-----------------------recover stresses at nip integrating points--------
nip=1; DEALLOCATE(points,weights)
ALLOCATE(points(nip,ndim),weights(nip))
CALL sample(element,points,weights)
!
WRITE(11,'(A,I2,A)')" The integration point (nip=",nip,") stresses are:"
!
WRITE(11,'(A,A)')" Element x-coord
y-coord",
&
!
"
sig_x
sig_y
tau_xy"
elements_5: DO iel=1,nels
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(:ndof))
gauss_pts_2: DO i=1,nip
CALL shape_fun(fun,points,i); CALL shape_der(der,points,i)
gc=MATMUL(fun,coord); jac=MATMUL(der,coord); CALL invert(jac)
deriv=MATMUL(jac,der); CALL beemat(bee,deriv)
sigma=MATMUL(dee,MATMUL(bee,eld))
