eval(demosdt('echoon'));format short e
clear variables global; comgui('close all');

% Demonstration of the model updating capabilities of the SDT
% for the case of the GARTEUR SM-AG-19 Testbed.


% initialize by loading results of the GARTFE and GARTTE demos and building
% the correlation matrix as in the GARTCO demo

load sdt_gart FEnode FEelt mdof md0 f0 sdof IIpo IIres pl il
femesh('selgroup1:10');femesh plotel0; cf=feplot;fecom(';sub1;view3');

% Now build a parametrized model based on the upcom interface

if ~exist('FileName','var')|(exist(FileName)~=2 & exist([FileName '.mat'])~=2)
 FileName = tempname;
end
warning(sprintf(['The superelement will be stored ' ...
 '%s.mat\n REMEMBER to delete it\n'],FileName));

clear global Up;
upcom(['load ' FileName]);
femesh('selgroup1:10'); % don't take group 11 which is the experimental mesh
upcom('set nominal t 6',FEnode,FEel0,pl,il)

% renumber md0 for later comparison
ind=fe_c(Up.DOF,mdof,'ind');md0(ind,:)=md0;mdof = Up.DOF;
cf.def={md0,Up.DOF,[f0]};

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

upcom('parstack reset');
upcom('parstackadd k','Tail','group3');
upcom('parstackadd t','Constrained Layer','group6');
upcom info


% now compute the response assuming that the thickness of group 6 elements
% is 1cm instead of 1.1.

upcom('parcoef',[1 1/1.1]);
def2 = upcom('compute mode full 6 20 1e3');
cf.def(2)=def2;fecom(';ch7;show2def')

% Note how much this affects the frequency of many modes.

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

% To use the update algorithms you should define the sensor configuration
% matrix with the same DOFs as those of the 'Up' superelement. Thus following
% the steps described in GARTTE

model=struct('Node',FEnode,'Elt',FEelt,'DOF',mdof);sens=fe_sens('model',model);
sens.tdof=sdof; sens=fe_sens('arigid',sens);

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

% compute reduced basis minimum residual expansion

upcom('parcoef',[1 1]); [m,k]=upcom('assemble');
st = sprintf('dynamic %g',f0(7)/2);
T = [fe_reduc(st,m,k,mdof,(sens.cta)') md0];
T = fe_norm(T,m,k,0); 
mdex3 = fe_exp(IIres.',IIpo(:,1)*2*pi,sens,m,k,mdof,T,'mdre');
feplot('initdefb',mdex3,Up.DOF,IIpo(:,1));

% compute and display strain energy distribution for the expanded modeshapes

Eres = upcom('eners',mdex3);
% energy levels in percent of total
Eres.data = Eres.data*diag(sum(Eres.data).^(-1));  
fecom('ch1');feplot('ColorDataElt',Eres);
colormap(jet(8));
fecom(';showpatch;colorbar');iimouse('ColorBarVHalf');

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

% compute dynamic residual and display the associated strain energies

res = ofact(k+1e3*m)\ (k*mdex3-m*mdex3*diag((IIpo(:,1)*2*pi).^2));
feplot('ColorDataElt',upcom('ener k',res));
fecom('sub');fecom('colorbar');

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

% Now do something more intelligent for parameter selection. Allow
% modifications of constraining layer thickness (group 6) and tail masses

upcom('parstack reset');
upcom('parstackadd k','Tail','group3');
upcom('parstackadd t','Constrained Layer','group6');
upcom('parstackadd m','Tail mass','group10');
upcom('optshift 1e3');upcom('optunit11');

% assuming that the parameter range defined in above is acceptable (otherwise 
% run a new series of 'ParStack' commands), the update based on frequencies is 
% simply obtained by

i2=1:8;[coef,md1,f1] = up_freq('basic',IIpo(i2,:),IIres(i2,:).',sens);

aa=[IIpo(:,1) f0(7:20) f1(7:20)];aa=[aa aa(:,2)./aa(:,1)-1 aa(:,3)./aa(:,1)-1];
fprintf(1,'tst %7.2f  ini%7.2f -> fin%7.2f (%4.1f -> %4.1f)\n',(aa*diag([1 1 1 100 100]))')

% Note that this result is not great because the parameter choice lowers
% modes 1-2 very much to improve the result for modes 3-5. Again, the
% up_freq algorithm should be taken as an example rather than a truly 
% efficient methodology.

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

% For large order models, you may want to use model reduction, for example
% you could use a basis combining nominal modes and sensitivities

upcom('parcoef',[1 1 1]);
[fsen,mdsen,mode,freq] = upcom('sens mode full',eye(3),7:20);

[m,k]=upcom('assemble');T = fe_norm([mdsen mode],m,k);
upcom('par red',[T])

upcom('parcoef',[1.1 1.2 1.26]);
t   =cputime;def1 = upcom('compute mode full 6 20 1e3');t(1)=cputime-t;
t(2)=cputime;def2 = upcom('compute mode reduced');t(2)=cputime-t(2);
cf.def(1)=def1; cf.def(2)=def2; fecom('ch7');
fecom('show2def');

fprintf(1,'Time for full order model (%.2f s) reduced (%.2f s)\n',t);

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

% This shows you how to compute FRFs efficiently (which is something you would
% in an FRF based update algorithm like UP_IXF).

IIw=linspace(1,60,1024)';b = fe_c(Up.DOF,58.03)';
upcom('parcoef',coef);
[IIxf,mode,freq] = upcom('compute modal full',0.01,b,b',IIw);
upcom('parcoef',[1 1 1]);
[IIxe,mode,freq] = upcom('compute modal full',0.01,b,b',IIw);
iicom(';submagpha;iixeon');

demosdt(['Input Do you want to delete ' FileName '.mat ?'], ...
   'Delete temp file ?','Yes','No','Yes', ...
   ['delete(''' FileName '.mat'');']);

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


%	Etienne Balmes  07/14/93, 04/25/00
%       Copyright (c) 1990-2003 by SDTools, All Rights Reserved.
%       $Revision: 1.6 $  $Date: 2005/06/06 13:52:12 $