%
disp(' ');
disp(' circular_plate_base_accel_fea.m  ver 1.8  October 8, 2012 ');
disp(' ');
disp(' Tom Irvine   Email: tom@vibrationdata.com');
disp(' ');
%
disp(' This program calculates the bending mode shapes and frequencies ');
disp(' of a circular plate using the Specht triangular element. ');
disp(' ');
disp(' It also calculates the frequency response functions ');
disp(' for base excitation.      ');
disp(' ');
disp(' A reference finite element model is used. ');
disp(' ');
disp(' Reference:  B. Specht, Modified Shape Functions for the Three-Node');
disp(' Plate Bending Element Passing the Patch Test, International'); 
disp(' Journal for Numerical Methods in Engineering, 1988. ');
disp(' ');
%
close all;
%
clear ModeShapes;
clear fn;
clear mass;
clear stiff;
clear m;
clear mass_local;
clear stiff_local;
clear mglobal_unc;
clear kglobal_unc;
clear A;
clear cbar;
clear b1;
clear b2;
clear b3;
clear c1;
clear c2;
clear c3;
clear fn;
clear omega;
clear E;
clear rho;
clear v;
clear D;
clear THM_elem;
%
fig_num=1;
%
[term_integ]=triangular_plate_specht_integrated_terms();
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Assume uniform plate materials and thickness
%
iu=1;
%
disp(' ');
disp(' Enter the diameter (inch) ');
diam=input(' ');
%
[amass,tmass,ins,imeth,thick,D,rho,v,h]=plate_materials_mass(iu);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear nodex;
clear nodey;
clear node1;
clear node2;
clear node3;
clear node4;
%
%%   The node coordinate matrix must have two columns:  x(inch)  y(inch) 
%
[THM]=circular_plate_grid_reference(1);
%
THM=THM*diam/128;
%
sz=size(THM);
dof=3*sz(1);
number_nodes=sz(1);
%
nodex=THM(:,1);
nodey=THM(:,2);
%
hold off;
figure(fig_num);
fig_num=fig_num+1;
plot(nodex,nodey,'.');
for i=1:number_nodes
   string=num2str(i,'%d\n');
%   text(nodex(i),nodey(i),string);
end
xlabel('X');
ylabel('Y');
grid on;
%
max_ux=max(nodex);
max_uy=max(nodey);
%
min_ux=min(nodex);
min_uy=min(nodey);
%
ux=max_ux-min_ux;
uy=max_uy-min_uy;
%
avex=(max_ux+min_ux)/2;
avey=(max_uy+min_uy)/2;
%
ddd=[ux,uy];
dmax=max(ddd);
dmax=dmax*1.2;
%
xmin=avex-dmax/2;
xmax=avex+dmax/2;
ymin=avey-dmax/2;
ymax=avey+dmax/2;
%
axis([xmin,xmax,ymin,ymax]);
hold on;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
%%  The element connectivity matrix must have three columns 
%%    node 1   node 2   node 3   
%
%
[THM]=circular_plate_elements_reference(1);
%
THM_elem=THM;
sz=size(THM);
nelem=sz(1);
%
node1=THM(:,1);
node2=THM(:,2);
node3=THM(:,3);
%
%  coincident node check
%
tol=1.0e-08;
%
for i=1:nelem
    iflag=0;
    if( abs( nodex(node1(i))-nodex(node2(i)) ) < tol && abs( nodey(node1(i))-nodey(node2(i)) ) < tol )
        iflag=1;
    end
    if( abs( nodex(node1(i))-nodex(node3(i)) ) < tol && abs( nodey(node1(i))-nodey(node3(i)) ) < tol )
        iflag=1;
    end   
    if( abs( nodex(node2(i))-nodex(node3(i)) ) < tol && abs( nodey(node2(i))-nodey(node3(i)) ) < tol )
        iflag=1;
    end    
    if(iflag==1)
        out1=sprintf(' Element error:  %d',i);
        disp(out1);
    end
end
%
for i=1:nelem
    x=[ nodex(node1(i)) nodex(node2(i)) nodex(node3(i)) nodex(node1(i)) ];
    y=[ nodey(node1(i)) nodey(node2(i)) nodey(node3(i)) nodey(node1(i)) ];
    plot(x,y);
    hold on;
end
hold off;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' Unconstrained model statistics ');
out1=sprintf('\n  nodes=%d  elements=%d   degrees-of-freedom=%d  \n',...
                  number_nodes,nelem,dof);
disp(out1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
mglobal_unc=zeros(dof,dof);
kglobal_unc=zeros(dof,dof);
%
total_area=0;
%
progressbar;
%
for ijk=1:nelem
%
    progressbar(ijk/nelem);
%
    x1= nodex(node1(ijk));
    x2= nodex(node2(ijk));
    x3= nodex(node3(ijk));
    y1= nodey(node1(ijk));
    y2= nodey(node2(ijk));
    y3= nodey(node3(ijk));

    [stiff_local,mass_local,area]=triangular_plate_local(D,x1,x2,x3,y1,y2,y3,term_integ);
%
%
    n1=(3*node1(ijk))-2;
    n2=(3*node2(ijk))-2;
    n3=(3*node3(ijk))-2;
%
    pr=[n1 n1+1 n1+2  n2 n2+1 n2+2  n3 n3+1 n3+2 ];
    pc=pr;       
%
    for i=1:9
        for j=1:9
           mglobal_unc(pr(i),pc(j))= mglobal_unc(pr(i),pc(j))+mass_local(i,j);            
           kglobal_unc(pr(i),pc(j))= kglobal_unc(pr(i),pc(j))+stiff_local(i,j);
        end
    end
%
    total_area=total_area+area;
%
end
%
total_volume=total_area*h;
%
rho_s=rho;
if(ins==1)
%%%     disp(' 1= Add mass   ');
%%%     disp(' 2= Specify total mass'); 
    if(imeth==1)
%%%       amass=input(' Enter added mass(lbm) ');
          rho=rho_s+(amass/total_volume);
    else
%%%       tmass=input(' Enter total mass(lbm) ');
          rho=tmass/total_volume;
    end
end
total_mass=total_volume*rho;
%
mglobal_unc=mglobal_unc*rho*h;
%
%
progressbar(1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  BC
%
clear mass;
clear stiff;
clear mass_unc;
clear stiff_unc;
%
 mass=mglobal_unc;
stiff=kglobal_unc;
%
 mass_unc=mglobal_unc;
stiff_unc=kglobal_unc;    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
sz=size(mass);
dof=(sz(1));
num=dof;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
nx=nodex;
ny=nodey;
%
out1=sprintf('\n number of dofs =%d \n',num);
disp(out1);
%
[ea,nd,nem,idn]=circular_plate_drive_nodes(nx,ny);
%
[con_nodes,num_con_nodes,con_dof,num_con_dof]=...
                    circular_plate_constraint_nodes(nx,ny,nd,idn,nem); 
%
mass_unc=mass;
stiff_unc=stiff;
%
if(num_con_dof>0)
%
    con_dof=sort(con_dof,'descend');
%
    progressbar;
    for i=1:num_con_dof
        progressbar(i/num_con_dof);
        nn=con_dof(i);
%       
        mass(nn,:)=[];
        mass(:,nn)=[];
%
        stiff(nn,:)=[];
        stiff(:,nn)=[];
%
        for j=1:nem
            if(ea(j)>nn)
                ea(j)=ea(j)-1;
            end
        end
%
    end
    pause(0.5);
    progressbar(1);
%
end
%
%
num=max(size(mass)); % may have changed due to constraints
nff=num-nem;
%
[nff,ngw,Mwd,Mww,fn,omegan,ModeShapes,MST,part,T1,T2]=...
 partition_matrices_accel(num,nem,ea,mass,stiff);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select modal damping input method ');
disp('  1=uniform damping for all modes ');
disp('  2=damping vector ');
idm=input(' ');
%
disp(' ');
if(idm==1)
   disp(' Enter damping ratio ');
   dampr=input(' ');
   damp(1:num)=dampr;
else
   disp(' Modal Damping Vector ');
   damp = input(' Enter the vector:  ');
end
%
[acc_frf,fstart,fend,freq]=...
 partition_matrices_accel_frf(num,nff,nem,ngw,Mwd,Mww,fn,damp,omegan,...
                                          ModeShapes,MST,part,T1,T2);
%
disp(' circular_plate_base_fea_post ');
%                                            
[trans_acc_frf,trans_vel_frf,trans_disp_frf,trans_rel_disp_frf]=...
   circular_plate_base_fea_post(con_dof,num_con_dof,nodex,nodey,acc_frf,...
                                                 freq,fstart,fend,fig_num); 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
PF=part;
%
disp(' ');
disp(' Apply half-sine base input?  1=yes  2=no ');
isb=input(' ');
if(isb==1)
%
  [t,acc,rel_disp,abase]=...
    mdof_half_sine_core(fn,damp,PF,ModeShapes,T1,T2,ngw);  
%  
  [acc_half_sine,rel_disp_half_sine,fig_num]=...
      circular_plate_th_plot(t,acc,rel_disp,abase,...
                          nodex,nodey,con_dof,num_con_dof,fig_num);
end
%
disp(' ');
disp(' Apply arbitrary time history base input?  1=yes  2=no ');
iar=input(' ');
if(iar==1)
%  
  [t,acc,rel_vel,rel_disp,abase]=...
                mdof_arbit_base_core(fn,damp,PF,ModeShapes,T1,T2,ngw);
%  
  [acc_arb,rel_disp_arb,fig_num]=...
               circular_plate_th_plot(t,acc,rel_disp,abase,...
                          nodex,nodey,con_dof,num_con_dof,fig_num);    
end
%
%
disp(' ');
disp(' Output arrays: ');
disp(' ');
disp('   trans_acc_frf');
disp('   trans_vel_frf');
disp('   trans_disp_frf');
disp('   trans_rel_disp_frf');
%
if(isb==1)
    disp('   acc_half_sine');
    disp('   rel_disp_half_sine');       
end
%
if(iar==1)
    disp('   acc_arb');
    disp('   rel_disp_arb');       
end
%