%
disp(' ');
disp(' circular_plate_bending_fea.m  ver 1.6  September 12, 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(' 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);
%
disp(' ');
disp(' Select boundary condition ');
disp('  1=fixed along circumference ');  
disp('  2=simply-supported along circumference'); 
disp('  3=fixed at discrete points');
disp('  4=simply-supported at discrete points');
disp('  5=free ');
%
itype=input(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
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;    
%
if(itype ~= 5)
    [mass,stiff,pp,con]=circular_plate_bc_apply(mass,stiff,itype);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  eigen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(itype ~= 5)
    rigid=0;    
else
    rigid=1;    
end
%
[fn,omega,ModeShapes,MST,nv]=Generalized_Eigen_fea(stiff,mass,rigid);
%
mstore=mass;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(itype ~= 3)
%
    [ModeShapes]=mode_shape_correction(dof,pp,con,ModeShapes);
% 
end 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[~]=pf_emm_fea_post1(mass_unc,dof,fn,total_mass,ModeShapes);
%
[fig_num]=...
triangular_elements_plot_modes(nodex,nodey,fig_num,ModeShapes,THM_elem,fn,number_nodes);
%
clear sum;