Civil Engineering Reference
In-Depth Information
np_types,nr,nst=4,nxe,nye
REAL(iwp)::alpha,beta,bot,cg_tol,det,dlam,dsbar,dslam,d3=3.0_iwp,ff,
&
fftol,fnew,fstiff,lode_theta,ltol,one=1.0_iwp,ptot,sigm,
&
small_tol=1.e-5_iwp,tloads,tol,top,up,zero=0.0_iwp
CHARACTER(LEN=15)::element='quadrilateral'
LOGICAL::converged,cg_converged
!-----------------------dynamic arrays------------------------------------
INTEGER,ALLOCATABLE::nf(:,:),g(:),no(:),num(:),g_num(:,:),g_g(:,:),
&
etype(:)
REAL(iwp),ALLOCATABLE::acat(:,:),acatc(:,:),bee(:,:),bdylds(:),bload(:), &
caflow(:),coord(:,:),d(:),daatd(:,:),ddylds(:),dee(:,:),der(:,:), &
deriv(:,:),diag_precon(:),dl(:,:),dload(:),dsigma(:),eld(:),eload(:), &
elso(:),eps(:),g_coord(:,:),jac(:,:),km(:,:),loads(:),p(:),points(:,:),&
prop(:,:),qinc(:),qinva(:),qinvr(:),qmat(:,:),ress(:),rmat(:,:), &
sigma(:),storkm(:,:,:),stress(:),tensor(:,:,:),totd(:),u(:),val(:,:), &
vmfl(:),vmfla(:),vmflq(:),vmtemp(:,:),weights(:),x(:),xnew(:),
&
x_coords(:),y_coords(:)
!-----------------------input and initialisation--------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)nxe,nye,cg_tol,cg_limit,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),etype(nels),storkm(ndof,ndof,nels),&
coord(nod,ndim),stress(nst),dl(nip,nels),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),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))
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)
WRITE(11,'(A,I6,A)')"There are ",neq," equations"
ALLOCATE(loads(0:neq),bdylds(0:neq),totd(0:neq),ddylds(0:neq),p(0:neq), &
x(0:neq),xnew(0:neq),u(0:neq),diag_precon(0:neq),d(0:neq))
!---------------loop the elements to set up element data------------------
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
END DO elements_1; CALL mesh(g_coord,g_num,12)
CALL sample(element,points,weights); diag_precon=zero
!----starting element stiffness integration, storage and preconditioner---
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; storkm(:,:,iel)=km
DO k=1,ndof; diag_precon(g(k))=diag_precon(g(k))+km(k,k); END DO
END DO elements_2; diag_precon(1:)=one/diag_precon(1:)
!-----------------------read load weightings------------------------------
READ(10,*)loaded_nodes; ALLOCATE(no(loaded_nodes),val(loaded_nodes,ndim))
