elastic_rope_model.m

function model = elastic_rope_model(params)
%function model = elastic_rope_model(params)
%
% model according to the first example of the experiments in "A
% Reduced Basis Method for Paramterized Variational Inequalities"
% (B.Haasdonk, J.Salomon, B.Wohlmuth).
%
% params.xnumintervals: number of intervals, i.e. 200 for paper

% I.Maier, B. Haasdonk 30.05.2011

if nargin < 1
  params = [];
end;
if ~isfield(params,'xnumintervals');
  params.xnumintervals = 202;
end;
if ~isfield(params,'xrange');
  params.xrange = [0,1];
end;

model.xrange = params.xrange;
model.xnumintervals = params.xnumintervals;
model.decomp_mode = 0;
model.RB_generation_mode = 'enrichment_by_supremizers';
model.mu_names = {'mu_1','mu_2'};
model.mu_ranges = {[10,50]/200,[-0.05,0.5]};

model.set_mu = @set_mu_default;
model.get_mu = @get_mu_default;

% default mu-values
model.mu_1 = 21.7157/200;
model.mu_2 = 0.1111;
model.mu_fix = 30/200;

% function pointers
model.gen_model_data = @vi_gen_model_data;
model.detailed_simulation = @vi_detailed_simulation;
model.gen_detailed_data = @vi_gen_detailed_data;
model.gen_reduced_data = @vi_gen_reduced_data;
model.reduced_data_subset = @vi_reduced_data_subset;
model.rb_simulation = @vi_rb_simulation;
model.rb_reconstruction = @vi_rb_reconstruction;
model.plot_sim_data = @vi_plot_sim_data;
model.plot_detailed_data = @vi_plot_detailed_data;
model.get_rb_size = @(model,detailed_data) size(detailed_data.RB_U,2);
model.get_dofs_from_sim_data = @(sim_data) sim_data.U;
model.is_stationary = 1;
model.gridtype = 'onedgrid';

model.operators = @elastic_rope_operators;
model.get_inner_product_matrix = @standard_inner_product_matrix;
model.get_supremizer_matrix = @elastic_rope_supremizer_matrix;

% data functions
obstacle_coefficients = @(dummy,params) [1; params.mu_2];
obstacle_components = @(glob,params) ...
    {-0.2*(sin(pi*glob) - sin(3*pi*glob)) - 0.5, (glob-1/2)};
model.obstacle = @(glob,params) ...
    eval_affine_decomp_general(obstacle_components, ...
                               obstacle_coefficients,glob,params);
%%% inf-sup, coercivity and continuity constant
model.beta_b = 1;
model.alpha_a = @(model) min(model.mu_1,model.mu_fix);
model.gamma_a = @(model) max(model.mu_1,model.mu_fix);
% discretization parameter
%model.H = params.H;
model.enable_error_estimator = 0;

% parameter for the snaphots computation
model.RB_numintervals = [4, 4];

% parameter for making work with new descr
model.descr = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% auxiliary functions

function K = standard_inner_product_matrix(model,model_data)
%function K = standard_inner_product_matrix(model,model_data)
K = 1/model_data.dx * (-diag(ones(1,model.xnumintervals-2),-1) + ...
                       2*eye(model.xnumintervals-1) ...
                       -diag(ones(1,model.xnumintervals-2),1));

function [A,B,f,g] = elastic_rope_operators(model,model_data)
%function [A,B,f,g] = elastic_rope_operators(model,model_data)
%
% provides offline/online decomposition

val_fix = model.mu_fix;

if model.decomp_mode == 2
  A = [model.mu_1; val_fix];
  B = 1;
  f = 1;
  g = model.obstacle([0,0],model);
    
else
  
  dx = model_data.dx;
  H = model.xnumintervals -1;
  
  A_2 = (-diag(ones(1,H-1),-1) + 2*eye(H) - diag(ones(1,H-1),1)) / ...
        dx;
  A_1 = A_2;
  
  ind = find((model_data.grid.X >= 0.5) & (model_data.grid.X < (0.5+dx)));
  
  % if xnumintervals even => this is exactly the middle point

  if mod(model.xnumintervals,2) == 0
    
    A_1((size(A_1,1)*(ind-1)+ind+1):end) = 0;
    A_1(ind,ind) = (0.5-model_data.grid.X(ind-1)) /(dx^2);
    A_1(ind-1,ind-1) = (0.5-model_data.grid.X(ind-2)) /(dx^2);
    A_1(ind,ind-1) = -(0.5-model_data.grid.X(ind-1)) /(dx^2);
    A_1(ind-1,ind) = A_1(ind,ind-1);
    
    A_2(1:(size(A_2,1)*(ind-2)+ind-2)) = 0;
    A_2(ind,ind) = (model_data.grid.X(ind+1)-0.5) /(dx^2);
    A_2(ind-1,ind-1) = (model_data.grid.X(ind)-0.5) /(dx^2);
    A_2(ind,ind-1) = -(model_data.grid.X(ind)-0.5) /(dx^2);
    A_2(ind-1,ind) = A_2(ind,ind-1);
  else % if number of intervals odd : middle cell in both comp, * 0.5
    error('to be implemented for odd numbers!');
  end;
  
  if model.decomp_mode == 0
    A = model.mu_1 * A_1 + val_fix * A_2;
    f = -ones(H,1)*model_data.dx;
%    f = -ones(H,1);
    g = -model.obstacle(model_data.grid.X(2:end-1)',model);
  else %model.decomp_mode = 1
    A = {A_1, A_2};
    f = {-ones(H,1)*model_data.dx};
    gmin = model.obstacle(model_data.grid.X(2:end-1)',model);
    g = {-gmin{1},-gmin{2}};
  end;
  
  B = -eye(H);
  
end;