Civil Engineering Reference
In-Depth Information
CALL num_to_g(num,nf,g); g_num(:,iel)=num
g_coord(:,num)=TRANSPOSE(coord); g_g(:,iel)=g
IF(nband<bandwidth(g))nband=bandwidth(g)
END DO elements_1
WRITE(11,'(A,I5,A,/,A,I5,/,A,I5)')" There are",neq," equations",
&
" The half-bandwidth (including diagonal) is",nband+1
ALLOCATE(kb(neq,nband+1),mb(neq,nband+1),ua(0:neq),va(0:neq),
&
diag(0:neq),udiag(0:neq),w1(0:neq),y(0:neq,leig),v_store(0:neq,lalfa))
kb=zero; mb=zero; ua=zero; va=zero; eig=zero; jeig=0; x=zero; del=zero
nu=0; alfa=zero; beta=zero; diag=zero; udiag=zero; w1=zero; y=zero; z=zero
CALL sample(element,points,weights)
!----------------- element stiffness integration and assembly-------------
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; mm=zero
integrating_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)
call ecmat(ecm,fun,ndof,nodof)
mm=mm+ecm*det*weights(i)*prop(3,etype(iel))
END DO integrating_pts_1
CALL formkb(kb,km,g); CALL formkb(mb,mm,g)
END DO elements_2; CALL cholin(mb)
!------------------------------find eigenvalues---------------------------
DO iters=1,lalfa
call lancz1(neq,el,er,acc,leig,lx,lalfa,lp,iflag,ua,va,eig,jeig,neig,x,&
del,nu,alfa,beta,v_store)
IF(iflag==0)EXIT
IF(iflag>1)THEN; WRITE(11,'(A,I5)')
&
" Lancz1 is signalling failure, with iflag = ",
iflag; STOP
END IF
!--- iflag = 1 therefore form u + a * v ( candidate for ebe ) -----------
udiag=va; CALL chobk2(mb,udiag); CALL banmul(kb,udiag,diag)
CALL chobk1(mb,diag); ua=ua+diag
END DO
!--- iflag = 0 therefore write out the spectrum -------------------------
WRITE(11,'(A,I5,A/)')" It took",iters," iterations"
WRITE(11,'(3(A,E12.4))')" Eigenvalues in the range",el," and",er," are:"
WRITE(11,'(6E12.4)')eig(1:neig)
!------------------- calculate the eigenvectors --------------------------
IF(neig>10)neig=10
call lancz2(neq,lalfa,lp,eig,jeig,neig,alfa,beta,lz,jflag,y,w1,z,v_store)
!-------------------if jflag is zero calculate the eigenvectors ---------
IF(jflag==0)THEN
DO i=1,nmodes
udiag(:)=y(:,i); CALL chobk2(mb,udiag)
WRITE(11,'(" Eigenvector number",I4," is:")')i
WRITE(11,'(6E12.4)')udiag(1:)/MAXVAL(ABS(udiag(1:)))
END DO
ELSE
! lancz2 fails
WRITE(11,'(A,I5)')" Lancz2 is signalling failure with jflag = ",jflag
END IF
STOP
END PROGRAM p103
