It's neither a in-box application nor in-house code, conversely a public-house
.It is published in book and I want it to improve a little further.
Where the hell has gone these attachments this time. Something is ill-conditioned with this server. ok I upload them to different location.
Code is below:
Program 4.4 Analysis of elastic rigid-jointed frames using 2-node beam/rod elements in 2- or 3-dimensions.
USE main; IMPLICIT NONE
INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
INTEGER::fixed_freedoms,i,iel,k,loaded_nodes,ndim,ndof,nels,neq,nod=2,
nodof,nn,nprops,np_types,nr REAL(iwp):
enalty=l.0e2 0_iwp,zero=0.0_iwp
i dynami c arrays
INTEGER,ALLOCATABLE::etype{
,g
),g_g
,
,g_num
,
,kdiag
),nf
,
no
),node
),num{
,sense
) REAL(iwp),ALLOCATABLE::action
),coord
,
,eld{
,gamma
),g_coord
,:
km
,
,kv
),loads
) ,prop
,
,value
)
i input and initialisation
OPEN(10,FILE=»fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)nels,nn,ndim,nprops,np_types
IF(ndim==2)nodof=3; IF(ndim==3)nodof=6; ndof=nod*nodof
ALLOCATE(nf(nodof,nn),km(ndof,ndof),coord(nod,ndim),g_coord(ndim,nn),
eld(ndof),action(ndof),g_num(nod,nels),num(nod),g(ndof),gamma(nels),
g_g(ndof,nels),prop(nprops,np_types),etype(nels)) READ(10,*)prop; etype=l; IF(np_types>l)READ(10,*)etype; IF(ndim==3)READ(10,*)gamma READ (10 , *) g_coord; READ (10, *) g__num
nf=l; READ(10,*)nr,(k,nf
,k),i=l,nr); CALL formnf(nf); neq=MAXVAL(nf) ALLOCATE(kdiag(neq),loads(0:neq)); kdiag=0 i loop the elements to find global array sizes
STATIC EQUILIBRIUM OF STRUCTURES
129
elements_l: DO iel=l,nels
num=g__num
, iel) ; CALL num_to_g (num, nf, g) ; g_g
, iel) =g
CALL fkdiag(kdiag,g) END DO elements__l
D0 i=2,neq; kdiag(i)=kdiag(i)+kdiag(i-1); END DO; ALLOCATE(kv(kdiag(neq))) WRITE(11, ' (2(A,15)) ') &
« There are",neq," equations and the skyline storage is",kdiag(neq)
i global stiffness matrix assembly
jcv^zero
elements_2: DO iel=l,nels
num=g_num
, iel) ; coord=TRANSPOSE (g__coord
, num) )
CALL rigid_jointed(km,prop,gamma,etype,iel,coord); g=g_g
,iel)
CALL fsparv(kv,km,g,kdiag) END DO elements_2
i read loads and/or displacements
loads=zero; READ(10,*)loaded_nodes,(k,loads(nf
,k)),i=l,loaded_nodes)
READ(10,*)fixed_freedoms
IF(fixed_freedoms/=0)THEN
ALLOCATE (node (fixed__f reedoms) , no (f ixed_f reedoms) , &
sense (fixed_f reedoms) , value (fixed_f reedoms) )
READ(10,*)(node(i),sense(i),value(i),i=l,fixed_freedoms)
DO i=l,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
i equation solution
CALL sparin(kv,kdiag); CALL spabac(kv,loads,kdiag); loads(0)=zero WRITE(11, ' (/A) ') " Node Displacements and Rotation(s)" DO k=l,nn; WRITE(11,'(15,6E12.4)?)k,loads(nf
,k)); END DO
i retrieve element end actions
WRITE(11,'(/A)')" Element Actions" elements_3: DO iel=l,nels
num=g_num
,iel); coord=TRANSPOSE(g_coord
,num))
CALL rigid_jointed(km,prop,gamma,etype,iel,coord); g=g_g
,iel)
eld=loads(g); action=MATMUL(km,eld)
IF(ndim<3)THEN; WRITE(11,•(15,6E12.4)')iel,action; ELSE WRITE(11,« (I5,6E12.4)')iel, action(l: 6) WRITE(11,'(A,6E12.4)')» ",action(7:12)
END IF END DO elements_3 STOP END PROGRAM p44
