eval(demosdt('echoon'))

% This demonstration gives an example of component mode synthesis using
% superelements.
%
% This is a fairly advanced demonstration so take your time ...
%
% d_cms illustrates the use of direct methods to achieve the same
% result.  See the second section in doc('sdt/trd') for the HTML help.

% The following lines initialize the problem (see FEMESH or the associated
% demos for more details on the commands).


femesh('reset');
FEnode=[1 0 0 0  0 .015 0;2 0 0 0 0 .028 0;3 0 0 0  0 .032 0];FEelt=[];
femesh(';objectbeamline1 2;divide 4;rotatesel 0 10 0 0 1;rev 4 0 -20 0 0 1');
femesh('addsel');
femesh(';objectbeamline2 3;divide 2;rotatesel 0 60 0 0 1;rev 24 0 -120 0 0 1');
femesh('addsel');
i1=femesh('findnode x>-.003 & x<.003 & y>.03');
femesh('objectbeamline',i1([2 1 3]));femesh('extrude 10 0 .002 0');
femesh('addsel');
femesh(';joingroup1 2 3;selgroup1')
femesh(';extrude 1 0.0 0.0 2e-3');FEelt=FEel0;
femesh(';orient;plotelt');fecom(';guione;sub;view3;showline');

fesuper('new blade',FEnode,FEelt);
global SEblade

SEblade=fe_case(SEblade,'fixdof','Fixed root','rad <.02 0');
SEblade.pl=m_elastic('dbval 1 Aluminum');
SEblade.il=p_solid('dbval 1 Full');
SEblade=fe_mk(SEblade)

SEblade.DOF=fe_case(SEblade,'gettdof'); % from now on ignore original model.DOF


% make the superelement generic, define reference orientation, create the two
% missing blades

FEel0=fesuper('make blade generic');
fesuper('set blade ref 1   1 2 1 10   1 2 1 10');
femesh(';addsel;rotatesel 0 120 0 0 1;addsel;');
femesh(';rotatesel 0 120 0 0 1;addsel;');
femesh('selgroup2:4');femesh plotelt;

cf=feplot;comgui('wplace tr');

[m,k,mdof]=fe_mk(FEnode,FEel0,[],[]);
cf.def=fe_eig({m,k,mdof},[6 20 0 11]);


demosdt('pause'); %----------------------------------------------------------

% now get model matrices for the first blade, find interface DOFs with other
% blades, and do a Craig-Bampton reduction (fixed interface normal modes and
% constraint modes). 

femesh('selgroup2');[m,k,mdof]=fe_mk(FEnode,FEel0,[],[]);
idof = fe_c(mdof,femesh('findnode group3:4 & group2'),'dof');

[T,sdof] = fe_reduc('CraigBampton 15 0 11',m,k,mdof,idof);

% project superelement and compute the modes of the 3 blade disk

fesuper('set blade tr',T,mdof,sdof);

femesh('selgroup2:4'); red_mod=femesh('model0');

[mr,kr,rdof]=fe_mk(red_mod);
[md2,f2]=fe_eig(mr,kr);

figure(1);comgui('wplace br');
plot([cf.def.data(:,1) f2(1:20)/2/pi]);legend('2016 DOF','99 DOF');
ylabel('Frequency');setlines([],{'-','--'},'+o');

% Clearly Component Mode Synthesis works very well for this example
% but feplot does not (yet) know how to reexpand to the full order model
% size. You could display the motion of the first blade

demosdt('pause'); %----------------------------------------------------------

[T,mdof,adof]=fesuper('get blade tr');
% make element DOF actual and not generic
ind=find(adof<0); eltid=1; adof(ind)=adof(ind)+1-eltid;  
fecom('sel1group1');
cf.def={T*md2(fe_c(rdof,adof,'ind'),:),mdof,f2};

% but the distinction between generic and actual element DOFs is difficult.

demosdt('pause'); %----------------------------------------------------------

% The other method is to define a visualization mesh (a tedious but useful
% approach). Here for example, you could do

i1 = [28 30 29 155 152 151 28 0 93:95 220 218 219 93 ...
   0 152 171 186 233 249 0 137 186 202 218];
femesh('object beamline',i1);fesuper('set blade patch',FEel0);
femesh('selgroup2:4');femesh('plotel0');

% Here we want to keep 15 fixed interface modes but associate them to
% the nodal DOFs of the 5 non-interface nodes used for the visualization mesh

fesuper('set blade dof',[idof;fe_c(mdof,[171 202 186 233 249]','dof')]);

idof =fesuper('get blade dof');T1 = T*fe_coor(fe_c(mdof,idof)*T,2);

fesuper('set blade k 1 1',T1'*k*T1); fesuper('set blade k 2 2',T1'*m*T1);
[mr,kr,rdof]=fe_mk(red_mod); [md3,f3]=fe_eig(mr,kr);

% The results with the two reduced bases are identical
% but the one defined on nodes is easier to display
r1=max(abs([f2(1:20)./f3(1:20)-1]));
fprintf('\nMaximum relative difference on frequencies %.5e\n',r1);

cf=feplot;cf.def={md3,rdof,f3};fecom('scd.01');


%------------------------------------------------------------------------
eval(demosdt('echooff'))

%       Etienne Balmes  07/28/96, 28/04/03
%       Copyright (c) 1996-2003 by SDTools
%       All Rights Reserved.
%       $Revision: 1.4 $  $Date: 2005/10/26 09:51:08 $