Civil Engineering Reference
In-Depth Information
CHARACTER(LEN=15)::element='quadrilateral'; LOGICAL::cg_converged
!-----------------------dynamic arrays------------------------------------
INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),nf(:,:),no(:),
&
num(:)
REAL(iwp),ALLOCATABLE::al(:),ans(:),bee(:,:),c(:,:),coord(:,:),d(:),
&
dee(:,:),der(:,:),derf(:,:),deriv(:,:),derivf(:,:),diag_precon(:),
&
eld(:),fun(:),funf(:),gc(:),g_coord(:,:),jac(:,:),kay(:,:),kd(:,:),
&
ke(:,:),km(:,:),kc(:,:),lf(:,:),loads(:),p(:),points(:,:),prop(:,:),
&
sigma(:),storkd(:,:,:),storke(:,:,:),u(:),val(:,:),vol(:),volf(:,:),
&
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,np_types
CALL mesh_size(element,nod,nels,nn,nxe,nye)
ALLOCATE(dee(nst,nst),points(nip,ndim),coord(nod,ndim),derivf(ndim,nodf),&
jac(ndim,ndim),kay(ndim,ndim),der(ndim,nod),deriv(ndim,nod), &
derf(ndim,nodf),funf(nodf),bee(nst,ndof),km(ndof,ndof),kc(nodf,nodf), &
g_g(ntot,nels),ke(ntot,ntot),kd(ntot,ntot),c(ndof,nodf),fun(nod), &
x_coords(nxe+1),y_coords(nye+1),vol(ndof),nf(nodof,nn),g(ntot), &
volf(ndof,nodf),g_coord(ndim,nn),g_num(nod,nels),num(nod),weights(nip),&
storke(ntot,ntot,nels),storkd(ntot,ntot,nels),etype(nels),
&
prop(nprops,np_types),gc(ndim),sigma(nst),eld(ndof))
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
READ(10,*)x_coords,y_coords; READ(10,*)dtim,nstep,theta,npri,nres
nf=1; READ(10,*)nr,(k,nf(:,k),i=1,nr); CALL formnf(nf); neq=MAXVAL(nf)
ALLOCATE(loads(0:neq),ans(0:neq),p(0:neq),x(0:neq),xnew(0:neq),u(0:neq), &
diag_precon(0:neq),d(0:neq))
READ(10,*)loaded_nodes; ALLOCATE(no(loaded_nodes),val(loaded_nodes,ndim))
READ(10,*)(no(i),val(i,:),i=1,loaded_nodes)
READ(10,*)nlfp; ALLOCATE(lf(2,nlfp))
READ(10,*)lf; nls=FLOOR(lf(1,nlfp)/dtim); IF(nstep>nls)nstep=nls
ALLOCATE(al(nstep)); CALL load_function(lf,dtim,al)
!-----------------------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,'x')
g(1:15:2)=nf(1,num(:)); g(2:16:2)=nf(2,num(:)); g(17:)=nf(3,num(1:7:2))
g_num(:,iel)=num; g_coord(:,num)=TRANSPOSE(coord); g_g(:,iel)=g
END DO elements_1; CALL mesh(g_coord,g_num,12)
WRITE(11,'(A,I5,A)')" There are",neq," equations"
CALL sample(element,points,weights); diag_precon=zero
!----------element matrix integration, storage and preconditioner---------
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)
