eval(demosdt('echoon'))

% THIS IS A DEMO FOR ADVANCED USERS WHO UNDERSTAND CMS
%
% This demonstration gives an example of component mode synthesis using direct
% methods. d_cms2 illustrates the use of superelements do achieve the same
% result (this is really a more general way of doing it).
%
% See the second section in doc('sdt/trd') for the HTML help associated to
% this demo.

% 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  1 0 0;2 0 0 0 .5 0 0;3 0 0 0 .0 0 0;4 0 0 0 .0 0 .1; ...
        6 0 0 0  -.5 0 0;7 0 0 0 -1.0 0 0;8 0 0 0 0 0 -.1];
FEelt=[];
femesh(';objectbeamline 1:4;extrude 3 0 .5 0;addsel;');
femesh(';objectbeamline 8 3 6 7;extrude 3 0 .5 0;addsel;');
femesh('plotelt');fecom(';view(3);colordatagroup;show patch');
model=femesh;
model.pl = [1 fe_mat('m_elastic','SI',1) 190e9 .30 7800  (190e9/2/(1+.29))];
model.il=[1 fe_mat('p_shell','SI',1) 0 0 0 .01];
femesh(model);

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

% - define the max frequency model to 50 Hz
% - find the interface nodes and the DOFs used for a model containing
%   two components

IntNode = femesh('find node group1 & group2');
femesh('selgroup1');mo1=model;mo1.Elt=FEel0;
femesh('selgroup2');mo2=model;mo2.Elt=FEel0;

% - reduce component 1 using the Craig-Bampton Basis (constraint + 
%   fixed-interface modes)

% manual procedure
 mo1=fe_case(mo1,'fixdof','Fixed Interface',IntNode);
 d1=fe_eig(mo1,[5 50.1 1e3]);
 %Assemble model without boundary conditions
 [m1,k1,mdof]=fe_mknl(mo1,'NoT');
 tc1=fe_reduc('static',m1,k1,mdof,IntNode);
 T1 = struct('def',[tc1 d1.def],'DOF',mdof);

% Direct approach
 mo1=fe_case(mo1,'reset','DofSet','Interface',IntNode);
 RED=fe_reduc('craigbampton 50.1',mo1)  % cut-off frequency at 50.1 Hz
 % RED.TR.def and T1 are directly comparable

% - reduce component 2 using Rubin's basis (attachment + free-interface modes)

 %Assemble model without boundary conditions
 [m2,k2,mdof]=fe_mknl(mo2,'NoT');
 d2=fe_eig({m2,k2,mdof},[6 20 1e3]);
 ta2=fe_reduc('flex',m2,k2,mdof,IntNode,fe_c(mdof,IntNode(1),'dof'));
 T2 = struct('def',[ta2 d2.def(:,find(d2.data(:,1)<50*2*pi))],'DOF',mdof);

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

% - define the output shape matrix associated to the interface nodes

model.DOF=feutil('getdof',model); % combined model DOFs
cint = fe_c(model.DOF,IntNode);
T1=feutilb('placeindof',model.DOF,T1);
T2=feutilb('placeindof',model.DOF,T2);

% - use fe_coor to find a basis of the kernel the constraint that relative
%   motion should be equal to zero

T = [T1.def (T2.def-cint'*(cint*T2.def))]* ...
   fe_coor(cint*[T1.def -T2.def]);

% - do the coupled prediction by projecting the full order model on the
%   kernel

[m12,k12,mdof]=fe_mknl(model,'NoT');
[mdr,fr]=fe_norm(T,m12,k12);
def_cms=struct('def',mdr,'DOF',model.DOF,'data',fr/2/pi);

% which gives the same result than the longer
% [mdr,fr]=fe_eig(T'*(m1+m2)*T,T'*(k1+k2)*T,[2 20 1e3]);mdr = T*mdr;

% compute exact response and compare

def_full=fe_eig({m12,k12,model.DOF},[6 20 1e3]);
cf=feplot; cf.def(1)=def_cms; cf.def(2)=def_full;
fecom(';showline;show2def;scalematch; ch7');
ii_mac(def_full,def_cms,'inda',7:20,'indb',7:20,'mac errorplot -fig3')

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

% The following repeats the problem solved above but starts by duplicating
% interface nodes to that the element DOFs are disjoint

IntNode=femesh('unjoin 1 2');

model=femesh;model.DOF=feutil('getdof',model);
[m12,k12,mdof]=fe_mknl(model,'NoT');
femesh('selgroup1');mo1=model;mo1.Elt=FEel0;
femesh('selgroup2');mo2=model;mo2.Elt=FEel0;

 mo1=fe_case(mo1,'reset','DofSet','Interface',IntNode(:,1));
RED1=fe_reduc('craigbampton 3',mo1)  % three modes

 mo2=fe_case(mo2,'reset','DofSet','Interface',IntNode(:,2));
RED2=fe_reduc('craigbampton 3',mo2)  % tree modes

T1=feutilb('placeindof',model.DOF,RED1.TR);
T2=feutilb('placeindof',model.DOF,RED2.TR);

cint1 = fe_c(model.DOF,IntNode(:,1));
cint2 = fe_c(model.DOF,IntNode(:,2));
T = [T1.def T2.def]*fe_coor([cint1*T1.def -cint2*T2.def]);
[mdr,fr]=fe_norm(T,m12,k12);
def_cms=struct('def',mdr,'DOF',model.DOF,'data',fr/2/pi);

ind=7:20;figure(5);plot([def_full.data(ind) def_cms.data(ind) ],':+');
set(gca,'ylim',[4 max(def_full.data(ind))]);
legend('True','CMS');xlabel('Mode #');ylabel('Freq');setlines

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

%	Etienne Balmes  
%       Copyright (c) 1990-2006 by SDTools, All Rights Reserved
%       $Revision: 1.9 $  $Date: 2006/03/22 17:24:09 $