Civil Engineering Reference
In-Depth Information
!-----------------------dynamic arrays------------------------------------
INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),nf(:,:),no(:),
&
node(:),num(:),sense(:)
REAL(iwp),ALLOCATABLE::bee(:,:),coord(:,:),d(:),dee(:,:),der(:,:), &
deriv(:,:),diag_precon(:),eld(:),fun(:),gc(:),g_coord(:,:),g_pmul(:,:),&
g_utemp(:,:),jac(:,:),km(:,:),loads(:),p(:),points(:,:),prop(:,:), &
sigma(:),store(:),storkm(:,:,:),u(:),value(:),weights(:),x(:),xnew(:), &
x_coords(:),y_coords(:),z_coords(:)
!-----------------------input and initialisation--------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
read(10,*)nod,nxe,nye,nze,nip,cg_tol,cg_limit,np_types
CALL mesh_size(element,nod,nels,nn,nxe,nye,nze); ndof=nod*nodof
ALLOCATE(nf(nodof,nn),points(nip,ndim),dee(nst,nst),coord(nod,ndim), &
jac(ndim,ndim),der(ndim,nod),deriv(ndim,nod),fun(nod),gc(ndim), &
bee(nst,ndof),km(ndof,ndof),eld(ndof),sigma(nst),g_coord(ndim,nn), &
g_num(nod,nels),weights(nip),num(nod),g_g(ndof,nels),x_coords(nxe+1), &
g(ndof),y_coords(nye+1),z_coords(nze+1),etype(nels),g_pmul(ndof,nels), &
prop(nprops,np_types),storkm(ndof,ndof,nels),g_utemp(ndof,nels))
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
READ(10,*)x_coords,y_coords,z_coords
nf=1; READ(10,*)nr,(k,nf(:,k),i=1,nr); CALL formnf(nf); neq=MAXVAL(nf)
WRITE(11,'(A,I5,A)')" There are",neq," equations"
ALLOCATE(p(0:neq),loads(0:neq),x(0:neq),xnew(0:neq),u(0:neq),
&
diag_precon(0:neq),d(0:neq))
CALL sample(element,points,weights); diag_precon=zero
!----------element stiffness integration, storage and preconditioner------
elements_1: DO iel=1,nels
CALL hexahedron_xz(iel,x_coords,y_coords,z_coords,coord,num)
CALL num_to_g(num,nf,g); g_num(:,iel)=num
g_coord(:,num)=TRANSPOSE(coord); g_g(:,iel)=g
CALL deemat(dee,prop(1,etype(iel)),prop(2,etype(iel)))
num=g_num(:,iel); g=g_g(:,iel); coord=TRANSPOSE(g_coord(:,num)); 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_1
!-----------------------invert the preconditioner and get starting loads--
loads=zero; READ(10,*)loaded_nodes,(k,loads(nf(:,k)),i=1,loaded_nodes)
READ(10,*)fixed_freedoms
IF(fixed_freedoms/=0)THEN
ALLOCATE(node(fixed_freedoms),sense(fixed_freedoms),
&
value(fixed_freedoms),no(fixed_freedoms),store(fixed_freedoms))
READ(10,*)(node(i),sense(i),value(i),i=1,fixed_freedoms)
DO i=1,fixed_freedoms; no(i)=nf(sense(i),node(i)); END DO
diag_precon(no)=diag_precon(no)+penalty; loads(no)=diag_precon(no)*value
store=diag_precon(no)
END IF
diag_precon(1:)=one/diag_precon(1:); diag_precon(0)=zero
d=diag_precon*loads; p=d; x=zero; cg_iters=0
!-----------------------pcg equation solution-----------------------------
pcg: DO
cg_iters=cg_iters+1; u=zero
elements_2: DO iel=1,nels; g_pmul(:,iel)=p(g_g(:,iel)); END DO elements_2
