Civil Engineering Reference
In-Depth Information
reason, later programs in this Chapter, Programs 5.5 and 5.6, introduce solution methods
in which assembly of the global stiffness matrix is avoided entirely.
Program 5.4 General two- (plane strain) or three-dimensional analysis of elastic
solids.
PROGRAM p54
!-------------------------------------------------------------------------
! Program 5.4 General two- (plane strain) or three-dimensional analysis
! of elastic solids (optional gravity loading).
!-------------------------------------------------------------------------
USE main; USE geom; IMPLICIT NONE
INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
INTEGER::fixed_freedoms,i,iel,k,loaded_nodes,ndim,ndof,nels,neq,nip,nn, &
nod,nodof,nprops=3,np_types,nr,nst
REAL(iwp)::det,penalty=1.0e20_iwp,zero=0.0_iwp
CHARACTER(len=15)::element
!-----------------------dynamic arrays------------------------------------
INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),kdiag(:),nf(:,:), &
no(:),node(:),num(:),sense(:)
REAL(iwp),ALLOCATABLE::bee(:,:),coord(:,:),dee(:,:),der(:,:),deriv(:,:), &
eld(:),fun(:),gc(:),gravlo(:),g_coord(:,:),jac(:,:),km(:,:),kv(:),
&
loads(:),points(:,:),prop(:,:),sigma(:),value(:),weights(:)
!-----------------------input and initialisation--------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)element,nod,nels,nn,nip,nodof,nst,ndim,np_types; ndof=nod*nodof
ALLOCATE(nf(nodof,nn),points(nip,ndim),dee(nst,nst),g_coord(ndim,nn), &
coord(nod,ndim),jac(ndim,ndim),weights(nip),num(nod),g_num(nod,nels), &
der(ndim,nod),deriv(ndim,nod),bee(nst,ndof),km(ndof,ndof),eld(ndof),
&
sigma(nst),g(ndof),g_g(ndof,nels),gc(ndim),fun(nod),etype(nels),
&
prop(nprops,np_types))
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
READ(10,*)g_coord; READ(10,*)g_num
IF(ndim==2)CALL mesh(g_coord,g_num,12)
nf=1; READ(10,*)nr,(k,nf(:,k),i=1,nr); CALL formnf(nf); neq=MAXVAL(nf)
ALLOCATE(kdiag(neq),loads(0:neq),gravlo(0:neq)); kdiag=0
!-----------------------loop the elements to find global arrays sizes-----
elements_1: DO iel=1,nels
num=g_num(:,iel); CALL num_to_g(num,nf,g); g_g(:,iel)=g
CALL fkdiag(kdiag,g)
END DO elements_1
DO i=2,neq; kdiag(i)=kdiag(i)+kdiag(i-1); END DO; ALLOCATE(kv(kdiag(neq)))
WRITE(11,'(2(A,I5))')
&
" There are",neq," equations and the skyline storage is",kdiag(neq)
!-----------------------element stiffness integration and assembly--------
CALL sample(element,points,weights); kv=zero; gravlo=zero
elements_2: DO iel=1,nels
CALL deemat(dee,prop(1,etype(iel)),prop(2,etype(iel))); num=g_num(:,iel)
coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); km=zero; eld=zero
int_pts_1: DO i=1,nip
CALL shape_fun(fun,points,i); 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)
eld(nodof:ndof:nodof)=eld(nodof:ndof:nodof)+fun(:)*det*weights(i)
