module defconst implicit none real(8),parameter::ggg=9.8D0 real(8),parameter::pi=3.14159265358979323846D0 end module defconst program f90_EGV_VFRM !***************************************************************** !2D Frame Natural Frequency Analysis !***************************************************************** use defconst implicit none integer::i,j,k,l,m,kk,ne,n,ia,ja,iii character::strcom*50 !Comment integer::NODT !Number of nodes integer::NELT !Number of elements integer::MATEL !Number of material sets integer::KOX !Number of restricted nodes in x-direction integer::KOY !Number of restricted nodes in y-direction integer::KOZ !Number of restricted nodes in rotation integer,allocatable::kakom(:,:) !Element connectivity integer,allocatable::matno(:) !Material set number real(8),allocatable::wm(:,:) !work for material properties real(8),allocatable::Em(:) !Elastic modulus of element real(8),allocatable::AA(:) !Section area of element real(8),allocatable::AI(:) !Moment of inertia of element real(8),allocatable::gamma(:) !Unit weight of element real(8),allocatable::x(:) !x-coordinate of node real(8),allocatable::y(:) !y-coordinate of node integer,allocatable::nokx(:) !Restricted node number in x-direction integer,allocatable::noky(:) !Restricted node number in y-direction integer,allocatable::nokz(:) !Restricted node number in rotation integer,allocatable::mindex(:) !Index for treatment of restricted nodes real(8),allocatable::ak(:,:) !Element stiffness matrix real(8),allocatable::am(:,:) !Element mass matrix real(8),allocatable::hen(:,:) !Transform matrix real(8),allocatable::tak(:,:) !Total stiffness matrix real(8),allocatable::tam(:,:) !Total mass matrix real(8),allocatable::ev(:) !Eigenvalue real(8),allocatable::vec(:,:) !Eigenvector real(8),allocatable::work(:) !work for eigenvector integer::negv !Number of eigenvalues for output integer::nt !Number of degrees of freedom of all nodes (NODT x nhen) integer::nod=2 !Number of nodes per element integer::nhen=3 !Number of degrees of freedom per node integer::mm !Number of degrees of freedom of reduced FE equation real(8)::fact0=1.0D-30 !0 judgement character::fnameR*50,fnameW*50,fnameA*50 !Input & output file name character :: linebuf*1000,dummy*200 character :: fmt1*200,fmt2*200,fmt3*200,fmt4*200 real(4) :: t1,t2 !Input from command line call getarg(1,dummy) !Number of eigenvalues for output call getarg(2,fnameR) !Input data file call getarg(3,fnameW) !Output data file read(dummy,*) negv !Convert character 'dummy' to integer 'negv' !To start the measurement of calculation time t1=DIFT() !*************************** !data input !*************************** open(11,file=fnameR,status='old') !/*--------------------------------------------------------------------------------------*/ !/*Comment*/ !/*--------------------------------------------------------------------------------------*/ read(11,'(a)') strcom !------------------------------------------------------------------------------------------ !Basic values for analysis !------------------------------------------------------------------------------------------ read(11,*) NODT,NELT,MATEL,KOX,KOY,KOZ !------------------------------------------------------------------------------------------ !To declare array size !------------------------------------------------------------------------------------------ nt=NODT*nhen allocate(kakom(1:NELT,1:nod)) allocate(matno(1:NELT)) allocate(wm(1:MATEL,1:4)) allocate(Em(1:NELT)) allocate(AA(1:NELT)) allocate(AI(1:NELT)) allocate(gamma(1:NELT)) allocate(x(1:NODT)) allocate(y(1:NODT)) if(1<=KOX)allocate(nokx(1:KOX)) if(1<=KOY)allocate(noky(1:KOY)) if(1<=KOZ)allocate(nokz(1:KOZ)) allocate(mindex(1:nt)) allocate(ak(1:nod*nhen,1:nod*nhen)) allocate(am(1:nod*nhen,1:nod*nhen)) allocate(hen(1:nod*nhen,1:nod*nhen)) allocate(tak(1:nt,1:nt)) allocate(tam(1:nt,1:nt)) allocate(ev(1:nt)) allocate(vec(1:nt,1:nt)) allocate(work(1:nt)) !------------------------------------------------------------------------------------------ !material properties (Em,AA,AI,gamma) !------------------------------------------------------------------------------------------ do i=1,MATEL read(11,*) wm(i,1),wm(i,2),wm(i,3),wm(i,4) end do !------------------------------------------------------------------------------------------ !Element-nodes relationship, material set number (node-1,node-2,matno) !------------------------------------------------------------------------------------------ do ne=1,NELT read(11,*) kakom(ne,1),kakom(ne,2),matno(ne) do i=1,MATEL if(i==matno(ne))then Em(ne) =wm(i,1) AA(ne) =wm(i,2) AI(ne) =wm(i,3) gamma(ne)=wm(i,4) end if end do end do deallocate(wm) !------------------------------------------------------------------------------------------ !(x,y) coordinates !------------------------------------------------------------------------------------------ do i=1,NODT read(11,*) x(i),y(i) end do !------------------------------------------------------------------------------------------ !Restricted node number !------------------------------------------------------------------------------------------ if(00)then mm=mm+1 mindex(mm)=i end if end do do i=1,mm ia=mindex(i) do j=1,mm ja=mindex(j) tak(i,j)=tak(ia,ja) tam(i,j)=tam(ia,ja) end do end do !----------------------------------------------------- !Calculation of eigenvalues !----------------------------------------------------- call GJACOBI(mm,tak,tam,ev,vec,nt) !Treatment of restricted parts do k=1,mm do i=1,nt work(i)=0.0D0 end do do i=1,mm ia=mindex(i) work(ia)=vec(i,k) end do do i=1,nt vec(i,k)=work(i) end do end do !***************************************************** !To print results !***************************************************** open(12, file=fnameW, access='append') write(12,'(a)') '*eigen value & eigen vector' !Number of degree of freedom kk=mm if(negvfmk)fmk=abs(a(i,j)) if(abs(b(i,j))>fmg)fmg=abs(b(i,j)) end do end do epsk=fmk*fact epsg=fmg*fact do i=1,nf do j=1,nf vec(i,j)=0.0 end do vec(i,i)=1.0 end do kmax=5*nf*nf do k=1,kmax ip=1 iq=2 amax=a(1,2)**2.0D0/(a(1,1)**2.0D0+a(2,2)**2.0D0) do i=1,nf-1 do j=i+1,nf dam=a(i,j)**2.0D0/(a(i,i)**2.0D0+a(j,j)**2.0D0) if(amax