Civil Engineering Reference
In-Depth Information
no(:),num(:)
REAL(iwp),ALLOCATABLE::acat(:,:),acatc(:,:),bee(:,:),bdylds(:),bload(:), &
caflow(:),coord(:,:),daatd(:,:),ddylds(:),dee(:,:),der(:,:),deriv(:,:),&
dl(:,:),dload(:),dsigma(:),eld(:),eload(:),elso(:),eps(:),g_coord(:,:),&
jac(:,:),km(:,:),kv(:),loads(:),points(:,:),prop(:,:),qinc(:),qinva(:),&
qinvr(:),qmat(:,:),ress(:),rmat(:,:),sigma(:),stress(:),tensor(:,:,:), &
totd(:),val(:,:),vmfl(:),vmfla(:),vmflq(:),vmtemp(:,:),weights(:),
&
x_coords(:),y_coords(:)
!-----------------------input and initialisation--------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)nxe,nye,fftol,ltol,np_types
CALL mesh_size(element,nod,nels,nn,nxe,nye)
ALLOCATE(nf(nodof,nn),points(nip,ndim),weights(nip),g_coord(ndim,nn), &
x_coords(nxe+1),y_coords(nye+1),num(nod),dee(nst,nst), &
tensor(nst,nip,nels),g_g(ndof,nels),coord(nod,ndim),jac(ndim,ndim), &
der(ndim,nod),deriv(ndim,nod),g_num(nod,nels),bee(nst,ndof), &
km(ndof,ndof),eld(ndof),eps(nst),sigma(nst),bload(ndof),eload(ndof), &
elso(nst),g(ndof),vmfl(nst),qinvr(nst),dl(nip,nels),stress(nst), &
dload(ndof),caflow(nst),dsigma(nst),ress(nst),rmat(nst,nst), &
acat(nst,nst),acatc(nst,nst),qmat(nst,nst),qinva(nst),daatd(nst,nst), &
vmflq(nst),vmfla(nst),vmtemp(1,nst),prop(nprops,np_types),etype(nels))
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
READ(10,*)x_coords,y_coords
nf=1; READ(10,*)nr,(k,nf(:,k),i=1,nr); CALL formnf(nf); neq=MAXVAL(nf)
ALLOCATE(kdiag(neq),loads(0:neq),bdylds(0:neq),totd(0:neq),ddylds(0:neq))
!-----------------------loop the elements to find global arrays sizes-----
kdiag=0
elements_1: DO iel=1,nels
CALL geom_rect(element,iel,x_coords,y_coords,coord,num,'y')
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
!--------------starting element stiffness integration and assembly--------
elements_2: 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); 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(no(loaded_nodes),val(loaded_nodes,ndim))
READ(10,*)(no(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"
totd=zero; tensor=zero; dl=zero; ptot=zero
load_increments: DO iy=1,incs
ptot=ptot+qinc(iy); bdylds=zero; loads=zero; iters=0
DO i=1,loaded_nodes; loads(nf(:,no(i)))=val(i,:)*qinc(iy); END DO
!-----------------------plastic iteration loop----------------------------
