%
disp(' ');
disp(' triangular_specht.m  ver 1.6  September 7, 2012 ');
disp(' ');
disp(' Tom Irvine   Email: tom@vibrationdata.com');
disp(' ');
%
disp(' This program calculates the bending mode shapes and frequencies ');
disp(' using the Specht triangular element. ');
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 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;
%
[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;
%
disp(' ');
disp(' The node coordinate matrix must have two columns:  x(inch)  y(inch) ');
disp(' ');
clear THM;
THM = input(' Enter the node coordinate matrix name:  ');
%
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(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' The element connectivity matrix must have three columns ');
disp('   node 1   node 2   node 3   ');
disp(' ');
clear THM;
THM = input(' Enter the element connectivity matrix name:  ');
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_mass=total_area*rho*h;
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
%
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;
%
disp('   ');
disp(' Apply boundary conditions?  1=yes  2=no ');
ibc=input(' ');
%
if(ibc==1)
    disp(' ');
    disp(' The constraint matrix must have three four columns:  node number   TZ RX RY ');
    disp(' A value of 1 indicates a constraint.  A value of 0 indicates free. ');
    disp(' ');
    clear THM;
    THM = input(' Enter constraint matrix name:  ');
%
    sz = size(THM);
    nc=sz(1);
%
    clear con;
    ij=1;
%
    progressbar;
    for i=1:nc
        progressbar(i/nc);
%        
        nodec=THM(i,1);
%        
        if(THM(i,2)==1)
           con(ij)=3*nodec-2;
           ij=ij+1;
        end
%        
        if(THM(i,3)==1)
           con(ij)=3*nodec-1;
           ij=ij+1;
        end
%     
        if(THM(i,4)==1)
           con(ij)=3*nodec;
           ij=ij+1;
        end       
%
    end
    progressbar(1);
%
    clear length;
    pp=length(con);
    con=con';
    con=sort(con,1,'descend');
    con;
%
    progressbar;
    for ij=1:pp
        progressbar(ij/pp);
%        
         k=con(ij);
         sz=size(mass);
         mass(:,k)=[];
         mass(k,:)=[];
        stiff(:,k)=[];
        stiff(k,:)=[];        
%
    end
    progressbar(1);
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  eigen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(ibc==2)
    isd=1;
else
    disp(' ');
    disp(' Does the system have rigid-body modes? ');
    disp(' 1=yes  2=no ');
    isd=input(' ');
end
%
rigid=0;
%
if(isd==1)
    rigid=1;
end
%
[fn,omega,ModeShapes,MST,nv]=Generalized_Eigen_fea(stiff,mass,rigid);
%
disp(' ');
disp(' Eigenvalue calculation time');
toc
%
mstore=mass;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(ibc==1)
    [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;