Static Equilibrium of Linear Elastic Solids - Programming the Finite Element Method - page 196

END DO int_pts_1
CALL fsparv(kv,km,g,kdiag); gravlo(g)=gravlo(g)-eld*prop(3,etype(iel))
END DO elements_2
loads=zero; READ(10,*)loaded_nodes,(k,loads(nf(:,k)),i=1,loaded_nodes)
loads=loads+gravlo; READ(10,*)fixed_freedoms
IF(fixed_freedoms/=0)THEN
ALLOCATE(node(fixed_freedoms),sense(fixed_freedoms),

&

value(fixed_freedoms),no(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
kv(kdiag(no))=kv(kdiag(no))+penalty; loads(no)=kv(kdiag(no))*value
END IF
!-----------------------equation solution---------------------------------
CALL sparin(kv,kdiag); CALL spabac(kv,loads,kdiag); loads(0)=zero
IF(ndim==3)THEN; WRITE(11,'(/A)')" Node

x-disp y-disp z-disp"

ELSE; WRITE(11,'(/A)')" Node x-disp y-disp"

END IF
DO k=1,nn; WRITE(11,'(I5,3E12.4)')k,loads(nf(:,k)); END DO
!-----------------------recover stresses at element Gauss-points----------
!nip=1; DEALLOCATE(points,weights); ALLOCATE(points(nip,ndim),weights(nip))
!CALL sample(element,points,weights)
WRITE(11,'(/A,I2,A)')" The integration point (nip=",nip,") stresses are:"
IF(ndim==3)THEN
WRITE(11,'(A,/,A)')"

Element x-coord y-coord z-coord", &

" sig_x sig_y sig_z tau_xy tau_yz tau_zx"

ELSE; WRITE(11,'(A,A)') " Element x-coord y-coord", &

" sig_x sig_y tau_xy"

END IF
elements_3: 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); eld=loads(g)
int_pts_2: DO i=1,nip
CALL shape_der(der,points,i); CALL shape_fun(fun,points,i)
gc=MATMUL(fun,coord); jac=MATMUL(der,coord); CALL invert(jac)
deriv=MATMUL(jac,der); CALL beemat(bee,deriv)
sigma=MATMUL(dee,MATMUL(bee,eld))
IF(ndim==3)THEN; WRITE(11,'(I8,4X,3E12.4)')iel,gc
WRITE(11,'(6E12.4)')sigma
ELSE; WRITE(11,'(I8,2E12.4,5X,3E12.4)')iel,gc,sigma
END IF
END DO int_pts_2
END DO elements_3
IF(ndim==2)THEN; CALL dismsh(loads,nf,0.05_iwp,g_coord,g_num,13)
CALL vecmsh(loads,nf,0.05_iwp,0.1_iwp,g_coord,g_num,14)
END IF
STOP
END PROGRAM p54

New dynamic real arrays:
gravlo loads generated by gravity

Perusal of Programs 5.1 and 5.3 in this chapter will show that they are essentially
identical. The shape function and derivative subroutines shape fun and shape der ,
and the [ B ] and [ D ] subroutines beemat and deemat can all generate the appropriate