Civil Engineering Reference
In-Depth Information
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; 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,i7))')
&
" There are",neq," equations and the skyline storage is",kdiag(neq)
CALL sample(element,points,weights); kv=zero
!-----------------------element stiffness integration and assembly--------
elements_2: DO iel=1,nels
ddt=d4*(one+prop(3,etype(iel)))/(d3*prop(2,etype(iel)))
IF(ddt<dt)dt=ddt; CALL deemat(dee,prop(2,etype(iel)),prop(3,etype(iel)))
num=g_num(:,iel); coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); km=zero
gauss_pts_1: 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)
km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
END DO gauss_pts_1; CALL fsparv(kv,km,g,kdiag)
END DO elements_2
!-----------------------read load weightings and factorise equations------
READ(10,*)loaded_nodes; ALLOCATE(node(loaded_nodes),val(loaded_nodes,ndim))
READ(10,*)(node(i),val(i,:),i=1,loaded_nodes); CALL sparin(kv,kdiag)
!-----------------------load increment loop-------------------------------
READ(10,*)tol,limit,incs; ALLOCATE(qinc(incs)); READ(10,*)qinc
WRITE(11,'(/A)')" step
load
disp
iters"
oldis=zero; totd=zero; tensor=zero; ptot=zero
load_incs: DO iy=1,incs
ptot=ptot+qinc(iy); iters=0; bdylds=zero; evpt=zero
!-----------------------plastic iteration loop----------------------------
its: DO
iters=iters+1; loads=zero
WRITE(*,'(A,F8.2,A,I4)')" load",ptot," iteration",iters
DO i=1,loaded_nodes; loads(nf(:,node(i)))=val(i,:)*qinc(iy); END DO
loads=loads+bdylds; CALL spabac(kv,loads,kdiag)
!-----------------------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
CALL deemat(dee,prop(2,etype(iel)),prop(3,etype(iel)))
num=g_num(:,iel); coord=TRANSPOSE(g_coord(:,num))
g=g_g(:,iel); eld=loads(g); bload=zero
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); stress=sigma+tensor(:,i,iel)
CALL invar(stress,sigm,dsbar,lode_theta)
!-----------------------check whether yield is violated-------------------
f=dsbar-SQRT(d3)*prop(1,etype(iel))
IF(converged.OR.iters==limit)THEN; devp=stress; ELSE
IF(f>=zero)THEN
dq1=zero; dq2=d3/two/dsbar; dq3=zero
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
