%
%    partition_matrices  ver 1.1  September 24, 2012
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Partition Mass & Stiffness Matrices
%
%%%%%%%%%%%%%%%%%%%%% Mdd Kdd %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function[fn,omega,ModeShapes,MST,Mwd,Kwd,Mww,Kww,TT,T1,T2,KT,MT]=...
    partition_matrices(mass,stiff,num,nfree,nem,free_dof,enforced_dof,enforced_string)
%
disp(' ')
enforced_string
disp(' ')
%
Mdd=mass;
Kdd=stiff;
for i=1:nfree
    row=free_dof(i);
    col=row;
%
    Mdd(row,:)=[];
    Mdd(:,col)=[];   
%
    Kdd(row,:)=[];
    Kdd(:,col)=[];   
%
end
%
%%%%%%%%%%%%%%%%%%%%% Mdf Kdf %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Mdf=mass;
Kdf=stiff;
for i=1:nfree
    row=free_dof(i);
%
    Mdf(row,:)=[];
    Kdf(row,:)=[]; 
%
end
%
for i=1:nem
    col=enforced_dof(i);
%
    Mdf(:,col)=[];
    Kdf(:,col)=[]; 
%
end
%
%%%%%%%%%%%%%%%%%%%%% Mfd Kfd %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Mfd=Mdf';
Kfd=Kdf';
%
%%%%%%%%%%%%%%%%%%%%%%% Mff Kff %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Mff=mass;
Kff=stiff;
for i=1:nem
    row=enforced_dof(i);
    col=row;
%
    Mff(row,:)=[];
    Mff(:,col)=[];   
%
    Kff(row,:)=[];
    Kff(:,col)=[];   
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
I=eye(nem,nem);
T2=eye(nfree,nfree);
%
if(strcmp(enforced_string, 'disp')==1)
   invMff=pinv(Mff);
   T1=-invMff*Mfd;
else
   invKff=pinv(Kff);
   T1=-invKff*Kfd;
end
%
TT=zeros(num,num);
%
TT(1:nem,1:nem)=I;
TT(nem+1:num,1:nem)=T1;
TT(nem+1:num,nem+1:num)=T2;
%
MP=zeros(num,num);
MP(1:nem,1:nem)=Mdd;
MP(1:nem,nem+1:num)=Mdf;
MP(nem+1:num,1:nem)=Mfd;
MP(nem+1:num,nem+1:num)=Mff;
%
KP=zeros(num,num);
KP(1:nem,1:nem)=Kdd;
KP(1:nem,nem+1:num)=Kdf;
KP(nem+1:num,1:nem)=Kfd;
KP(nem+1:num,nem+1:num)=Kff;
%
MT=TT'*MP*TT
%
KT=TT'*KP*TT
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% [fn,omega,ModeShapes,MST]=Generalized_Eigen(KT,MT,1);
%
disp(' ');

Mwd=MT(nem+1:num,1:nem);
Kwd=KT(nem+1:num,1:nem);
Mww=MT(nem+1:num,nem+1:num);
Kww=KT(nem+1:num,nem+1:num);
%
disp(' Calculating eigenvalues... ');
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen(Kww,Mww,1);