multiscale_thermalblock_model.m

function model = multiscale_thermalblock_model(params)
%
% A thermalblock model, the parameters of which are 
% coupled to give overall params.Qa instead of params.B1*params.B2 
% many parameters.
%
% The coupling can be prescribed by a vector params.parameter_mapping
% e.g. for B1=B2=3 and desired Qa=4: 
% params.parameter_mapping = [4 2 3 1 2 1 3 1 1];
%
% If this mapping is not given, a deterministically 
% random ordering is prescribed
%
% a field Qa is required in params to indicate the number of
% desired components in a

% B. Haasdonk 7.3.2012

model = thermalblock_model(params);
model.base_model = model;
model.base_Qa = length(model.mu_names);
model.Qa = params.Qa;
Qa = params.Qa;

% deterministic random assignment
if isfield(params,'parameter_mapping')
  model.parameter_mapping = params.parameter_mapping;
else
  defaultStream = RandStream.getDefaultStream;
  savedState = defaultStream.State;
  RandStream.setDefaultStream(RandStream('mt19937ar','seed',100000));
  model.parameter_mapping = floor(Qa * rand(1,model.base_Qa))+1;
  defaultStream.State = savedState;
  RandStream.setDefaultStream(defaultStream);
end;
parameter_first_ids = [];
for i = 1:Qa
  parameter_first_ids = [parameter_first_ids, ...
                    min(find(model.parameter_mapping==i))];
end;

% assign new model parts:
model.parameter_first_ids = parameter_first_ids;
model.mu_names = model.mu_names(parameter_first_ids);
model.mu_ranges = model.mu_ranges(parameter_first_ids);
model.set_mu = @my_set_mu;
model.get_mu = @my_get_mu;
model.mus = model.mus(model.parameter_first_ids);
%model.operators = @my_operators;
model.name = 'multiscale_thermalblock';

% default params for demo_rb_gui
model.axis_tight = 1;
model.yscale_uicontrols = 1.0;

%keyboard;

multiscale_diffusivity_tensor_coefficients = @(dummy,params) ...
    params.get_mu(params);
model.global_diffusivity_tensor = @(glob,params) ...
    eval_affine_decomp_general(...
        @multiscale_diffusivity_tensor_components, ...
        multiscale_diffusivity_tensor_coefficients, glob,params);
model.diffusivity_tensor = @(varargin) ...
      discrete_volume_values(model.global_diffusivity_tensor,varargin{:});

model = elliptic_discrete_model(model);

function mu = my_get_mu(model)
mu = model.mus;

function model = my_set_mu(model,mu)
base_mu = mu(model.parameter_mapping);
model.base_model = model.base_model.set_mu(model.base_model,base_mu);
model.mus = mu;
%model = set_mu_default(model,mu);

%function [A,b] = my_operators(model,glob)
%%model.base_model.decomp_mode = model.decomp_mode;
%%[A,b] = model.base_model.operators(model.base_model,glob);
%[A,b] = model.operators(model,glob);
%if model.decomp_mode == 0
%  % do nothing
%elseif model.decomp_mode == 1
%  % reduce number of A components
%  A_comp_new = cell(model.Qa,1);
%  for q = 1:model.Qa
%%    A_comp_new{q} = zeros(size(A{q})); % initialize with 0;
%    A_comp_new{q} = 0*A{q}; % initialize with 0, sparse if required
%  end;
%  for q = 1:model.base_Qa
%    q2 = model.parameter_mapping(q);
%    A_comp_new{q2} = A_comp_new{q2} + A{q};
%  end;  
%  A = A_comp_new;
%else 
%  % select correct Qa coefficients
%  A = A(model.parameter_first_ids);     
%end;

function res = multiscale_diffusivity_tensor_components(glob,params)
% diffusion tensor: each row four entries a11,a_21,a_12,a_22. 
% a11(x)=a22(x) = mu_i if x in block i, a12=a21 = 0. 
% for each point in glob find global block number
% xi: range 0 ... B1-1
xi = floor(glob(:,1)*params.B1);
i = find(xi>=params.B1);
if ~isempty(i)
  xi(i) = params.B1-1; 
end;
% xi: range 0 ... B1-1
yi = floor(glob(:,2)*params.B2);
i = find(yi>=params.B2);
if ~isempty(i);
  yi(i) = params.B2-1; 
end;
block_index = yi*params.B1+xi+1;
mu_index = params.parameter_mapping(block_index);
zeroblock = zeros(size(glob,1),4);
%res = cell(1,params.number_of_blocks);
res = cell(1,params.Qa);
%for q = 1:params.number_of_blocks;
for q = 1:params.Qa;
  block = zeroblock;
  i = find(mu_index==q);
  if ~isempty(i)
    block(i,1) = 1;
    block(i,4) = 1;
  end;
  res{q} = block;
end;