Civil Engineering Reference
In-Depth Information
END DO elements_1; CALL mesh(g_coord,g_num,12)
DO i=2,neq; kdiag(i)=kdiag(i)+kdiag(i-1); END DO; ALLOCATE(kv(kdiag(neq)))
WRITE(11,'(2(A,I5))')
&
" There are",neq," equations and the skyline storage is",kdiag(neq)
CALL sample(element,points,weights); loads=zero; kv=zero
!-----------------------global matrix assembly----------------------------
elements_2: DO iel=1,nels
kay=zero; DO i=1,ndim; kay(i,i)=prop(i,etype(iel)); END DO
CALL deemat(dee,prop(3,etype(iel)),prop(4,etype(iel))); num=g_num(:,iel)
coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); km=zero; c=zero; kc=zero
gauss_points_1: DO i=1,nip
!-----------------------elastic solid contribution------------------------
CALL shape_der(der,points,i)
jac=MATMUL(der,coord); det=determinant(jac); CALL invert(jac)
deriv=MATMUL(jac,der); CALL beemat(bee,deriv)
vol(:)=bee(1,:)+bee(2,:)
km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
!-----------------------fluid contribution--------------------------------
CALL shape_fun(funf,points,i); CALL shape_der(derf,points,i)
derivf=MATMUL(jac,derf)
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; store_kc(:,:,iel)=kc
CALL formke(km,kc,c,ke,theta); CALL fsparv(kv,ke,g,kdiag)
END DO elements_2
!-----------------------factorise equations-------------------------------
CALL sparin_gauss(kv,kdiag)
!-----------------------time stepping loop--------------------------------
WRITE(11,'(/A,I5)')" Results at node",nres
WRITE(11,'(A)')
&
"
time
load
x-disp
y-disp
porepressure"
WRITE(11,'(5E12.4)')0.0,0.0,loads(nf(:,nres))
time_steps: DO j=1,nstep
tot_load=SUM(al(1:j)); time=j*dtim; ans=zero
elements_3: DO iel=1,nels
g=g_g(:,iel); kc=store_kc(:,:,iel)
phi0=loads(g(ndof+1:)); phi1=MATMUL(kc,phi0)
ans(g(ndof+1:))=ans(g(ndof+1:))+phi1
END DO elements_3
!-----------------------apply loading increment---------------------------
DO i=1,loaded_nodes; ans(nf(1:2,no(i)))=val(i,:)*al(j); END DO
!-----------------------equation solution---------------------------------
CALL spabac_gauss(kv,ans,kdiag); loads=loads+ans; loads(0)=zero
IF(j/npri*npri==j)WRITE(11,'(5E12.4)')time,tot_load,loads(nf(:,nres))
!-----------------------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_4: 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)
