demo_vi.m

Go to the documentation of this file.

function demo_vi()
%function demo_vi()
% function testing the implementation of the RB-Method for
% Parametrized Variational Inequalities.
% In particular 1D FEM model of a hanging elastic rope with different 
% elasticities in subdomains.

% B. Haasdonk, I.Maier 31.05.2011

help demo_vi;

params.H = 201;

model = elastic_rope_model(params);
disp('model created');
model_data = gen_model_data(model);

fn = 'detailed_data_vi.mat';

if 1
  % detailed simulation and plot
  disp('begin detailed simulation');
  sim_data = detailed_simulation(model,model_data);
  disp('plotting ...');
  plot_params.plot_title = 'detailed solution & obstacle';
  p = plot_sim_data(model,model_data,sim_data,plot_params);
  
  %g = model.obstacle(model_data.X,model);
  %figure(p),
  %hold on;
  %plot(model_data.X,g,'c','Linewidth',4);
  %title('detailed solution & obstacle');
  %hold off;
  
  % rb generation
  disp('compute snapshots');
  detailed_data = gen_detailed_data(model,model_data);
  %save('detailed_data_vi.mat','model','detailed_data');
  
  reduced_data = gen_reduced_data(model,detailed_data);
  
  % rb simulation and plot
  disp('begin reduced simulation');
  rb_sim_data = rb_simulation(model,reduced_data);
  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
  disp('plotting ...');
  plot_params.plot_title = 'reduced solution & obstacle';
  p = plot_sim_data(model,model_data,rb_sim_data,plot_params);
  
  % error computation
  disp('compute errors and norms');
  K = model.get_inner_product_matrix(model,model_data);
  U_err = sqrt((rb_sim_data.U-sim_data.U)' * K * (rb_sim_data.U- ...
                                                  sim_data.U));
  L_R = K \ (rb_sim_data.L-sim_data.L);
  L_err = sqrt((rb_sim_data.L-sim_data.L)' * L_R);
  disp(['error (U-UN): ',num2str(U_err)]);
  disp(['error (L-LN): ',num2str(L_err)]);
  
  U_norm = sqrt(sim_data.U' * K * sim_data.U);
  L_R = K \ sim_data.L;
  L_norm = sqrt(sim_data.L' * L_R);
  disp(['norm U: ',num2str(U_norm)]);
  disp(['norm L: ',num2str(L_norm)]);
  if ~exist(fn);
    save(fn,'model','detailed_data');
  end;
  
end % if 0

% demo_rb_gui
if 1
  load(fn);
  model.mu_fix = 30/200;
  plot_params = [];
  plot_params.no_lines = 0;
  plot_params.axis_tight = 1;
  demo_rb_gui(model,detailed_data,[],plot_params);
end