Civil Engineering Reference
In-Depth Information
END DO; ALLOCATE(kv(kdiag_l(neq)),mv(kdiag_r(neq)))
WRITE(11,'(A,I5,A,/,2(A,I5),A)')
&
" There are",neq," equations."," Skyline storage is",kdiag_l(neq),
&
" to the left, and",kdiag_r(neq)," to the right."
c1=one/dtim/dtim/beta; kv=zero; mv=zero
!-----------------------global stiffness and mass matrix assembly---------
elements_2: DO iel=1,nels
num=g_num(:,iel); coord=TRANSPOSE(g_coord(:,num))
CALL stiff4(km,coord,prop(1,etype(iel)),prop(2,etype(iel)))
g=g_g(:,iel); area=zero; mm=zero
gauss_pts_1: DO i=1,nip
CALL shape_der(der,points,i)
jac=MATMUL(der,coord); det=determinant(jac); area=area+det*weights(i)
CALL shape_fun(fun,points,i)
IF(mtype(iel)=='c')THEN; CALL ecmat(ecm,fun,ndof,nodof)
mm=mm+ecm*det*weights(i)*c1*prop(3,etype(iel))
END IF
END DO gauss_pts_1
IF(mtype(iel)=='l')THEN
CALL elmat(area,c1*prop(3,etype(iel)),mm)
CALL fsparv(kv,mm,g,kdiag_l); CALL fsparv(mv,mm-km,g,kdiag_r); ELSE
CALL fsparv(kv,km+mm,g,kdiag_l); CALL fsparv(mv,mm,g,kdiag_r)
END IF
END DO elements_2
!-----------------------initial conditions and factorise equations--------
x0=zero; d1x0=one; d2x0=zero; CALL sparin(kv,kdiag_l); time=zero
!-----------------------time stepping loop--------------------------------
WRITE(11,'(/A,I5))')" Result at node",nres
WRITE(11,'(A)')"
time
x-disp
y-disp"
WRITE(11,'(4E12.4)')time,x0(nf(:,nres))
timesteps: DO j=1,nstep
time=time+dtim; d1x1=x0+d1x0*dtim+d2x0*pt5*dtim*dtim*(one-two*beta)
CALL linmul_sky(mv,d1x1,x1,kdiag_r); CALL spabac(kv,x1,kdiag_l)
d2x1=(x1-d1x1)/dtim/dtim/beta
d1x1=d1x0+d2x0*dtim*(one-gamma)+d2x1*dtim*gamma
x0=x1; d1x0=d1x1; d2x0=d2x1
IF(j/npri*npri==j)WRITE(11,'(4E12.4)')time,x0(nf(:,nres))
END DO timesteps
STOP
END PROGRAM p115
New dynamic integer arrays:
kdiag l
diagonal term locations on the left
kdiag r
diagonal term locations on the right
New dynamic character array:
mtype
set to
'l'
for lumped mass and
'c'
for consistent
Programs 11.3 and 11.4 used implicit time-marching algorithms, and Program 11.7
described later in this chapter uses an explicit approach. The program described now com-
bines the methods of implicit and explicit time integration in a single program. Although not
“mesh free” this procedure should allow economical bandwidths of the assembled matrices.
The idea is (e.g. Key, 1980) that a mesh may contain only a few elements which have a
very small explicit stability limit and are therefore best integrated implicitly. The remainder
of the mesh can be successfully integrated explicitly at reasonable time steps.
