Civil Engineering Reference
In-Depth Information
until it assumes the shape of the free surface at convergence. Program 7.4 is a general
program for the solution of Laplace's equation over two- or three-dimensional domains.
The final Program 7.5 describes an element-by-element version of Program 7.4 avoiding
the need for global matrix assembly. A parallel version of this program is also described
in Chapter 12.
As the problems in this chapter involve scalar fields with just one unknown at each
node, the programming is simplified in that
nod
is always equal to
ndof
, so the latter
is not required. A further change from the solid mechanics applications of the preceding
chapters, is that the “zero freedoms” data previously introduced through the
nf
array has
been removed. In this chapter, all fixed boundary conditions, whether equal to zero or not,
are applied through the
fixed freedom
data. Since all nodal freedoms are therefore
retained in the assembly and analysis,
nn
is always equal to
neq
. In the interests of
consistency with other programs in the topic however,
neq
has been retained.
Program 7.1 One-dimensional analysis of steady seepage using 2-node line elements.
PROGRAM p71
!-------------------------------------------------------------------------
! Program 7.1 One dimensional analysis of steady seepage using
! 2-node line elements.
!-------------------------------------------------------------------------
USE main; IMPLICIT NONE
INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
INTEGER::fixed_freedoms,i,iel,k,loaded_nodes,nels,neq,nod=2,nn,nprops=1, &
np_types
REAL(iwp)::penalty=1.0e20_iwp,zero=0.0_iwp
!-----------------------dynamic arrays------------------------------------
INTEGER,ALLOCATABLE::etype(:),kdiag(:),g_num(:,:),node(:),num(:)
REAL(iwp),ALLOCATABLE::disps(:),ell(:),kp(:,:),kv(:),kvh(:),loads(:),
&
prop(:,:),value(:)
!-----------------------input and initialisation--------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)nels,nn,np_types; neq=nn
ALLOCATE(ell(nels),num(nod),prop(nprops,np_types),etype(nels),
&
kp(nod,nod),g_num(nod,nels),kdiag(neq),loads(0:neq),disps(0:neq))
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
READ(10,*)ell; READ(10,*)g_num; kdiag=0
!-----------------------loop the elements to find global arrays sizes-----
elements_1: DO iel=1,nels
num=g_num(:,iel); CALL fkdiag(kdiag,num)
END DO elements_1
DO i=2,neq; kdiag(i)=kdiag(i)+kdiag(i-1); END DO
WRITE(11,'(2(A,I5))')
&
" There are",neq," equations and the skyline storage is",kdiag(neq)
ALLOCATE(kv(kdiag(neq)),kvh(kdiag(neq))); kv=zero
!-----------------------global conductivity matrix assembly---------------
elements_2: DO iel=1,nels
CALL rod_km(kp,prop(1,etype(iel)),ell(iel))
num=g_num(:,iel); CALL fsparv(kv,kp,num,kdiag)
END DO elements_2; kvh=kv
!-----------------------specify boundary values---------------------------
loads=zero; READ(10,*)loaded_nodes,(k,loads(k),i=1,loaded_nodes)
READ(10,*)fixed_freedoms
IF(fixed_freedoms/=0)THEN
ALLOCATE(node(fixed_freedoms),value(fixed_freedoms))
