thermalblock_model_hp.m

Go to the documentation of this file.

function descr = thermalblock_model_hp(params)
% Thermal Block example similar as described in the book of 
% A.T. patera and G. Rozza (just one parameter more)
%
% i.e. 
% - div ( a(x) grad u(x)) = q(x)    on Omega
%                   u(x)) = g_D(x)  on Gamma_D
%     a(x) (grad u(x)) n) = g_N(x)  on Gamma_N
%                       s = l(u) linear output functional
%
% where Omega = [0,1]^2 is divided into B1 * B2 blocks 
% QA := B1*B2. The heat conductivities are given by mu_i:
%
% this is a modified thermalblock model fpr hp test
% with only 2 paramerts and 9
% paramert domains and
% changed range (.04:1) 
% also the location of mus is changed to:
%   -------------------------
%   | mu_1  | mu_2  | mu_1 |
%   -------------------------
%   | mu_2  | mu_1  | mu_2 |
%   -------------------------
%   | mu_1  | mu_2  | mu_1 |
%   -------------------------
%
% `Gamma_D =` upper edge
% `Gamma_N = boundary(Omega)\Gamma_D`
%
% `a(x) = mu_i(x) if x\in B_i`
% `q(x) = 0`
% `g_D  = 0` on top
% `g_N  = 1` on lower edge, 0 otherwise
% `l(u)` = average over lower edge

% S. Langhof and B. Haasdonk 26.2.2011
% I. Maier 26.04.2011

descr = lin_stat_model_default;
descr.name = 'thermalblock';
descr.dmodel_constructor = @ThermalBlock.DetailedModel;

descr.set_mu = @thermalblock_set_mu;
descr.get_mu = @(descr) descr.mus(:);

if nargin == 0
    params.B1=3;
    params.B2=3;
%    descr.default = '3x3 demo version';
end
% number of blocks in X and Y direction
descr.B1 = params.B1;
descr.B2 = params.B2;
descr.number_of_blocks = params.B1*params.B2;
descr.number_of_mu = 2;
% each block has a different heat conductivity
% upper left has conductivity 1
mu_names = {};
mu_ranges = {};
mu_range = [1/50,1];
mu_names = [mu_names,{['mu1']},{['mu2']}];%[mu_1,mu_2]

for p = 1:descr.number_of_mu
  
  mu_ranges = [mu_ranges,{mu_range}];
end;
descr.mu_names = mu_names;
descr.mu_ranges = mu_ranges;

%default values 1 everywhere
descr.mus = ones(descr.number_of_mu,1);


%%%%%%%%% set data functions
descr.has_diffusivity = 1;
descr.has_output_functional = 1;
descr.has_dirichlet_values = 1;
descr.has_neumann_values = 1;

% zero dirichlet values, i.e. 1 component, Q_dir = 1
dirichlet_values_coefficients = @(dummy,params) [0];
dirichlet_values_components = @(glob,params) {zeros(size(glob,1),1)};
descr.dirichlet_values = @(glob,params) ...
    eval_affine_decomp_general(dirichlet_values_components, ...
                               dirichlet_values_coefficients,glob,params);

% 1/0 neumann values depending on edge, non parametric, i.e Q_neu = 1;
neumann_values_coefficients = @(dummy,params) 1; % single component
descr.neumann_values = @(glob,params) ...
    eval_affine_decomp_general(@thermalblock_neumann_values_components, ...
                               neumann_values_coefficients, 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. 
diffusivity_tensor_coefficients = @(dummy,params) ...
 [params.mus(1),params.mus(2),params.mus(1),params.mus(2),...
   params.mus(1),params.mus(2),params.mus(1),params.mus(2),...
    params.mus(1)]; 
%params.mus(:);
descr.diffusivity_tensor = @(glob,params) ...
    eval_affine_decomp_general(...
        @thermalblock_diffusivity_tensor_components, ...
        diffusivity_tensor_coefficients, glob,params);

% only useful for detailed simulation or nonlinear outputs:
%descr.output_functional = @output_functional_boundary_integral;

descr.operators_output = @fem_operators_output_boundary_integral;
% output weight function: simply 1 on lower edge, 0 else => Ql = 1 component
output_function_components = @(glob,descr) {1.0*(glob(:,2)<eps)};
output_function_coefficients = @(glob,descr) 1;
descr.output_function = @(glob,params) ...
    eval_affine_decomp_general(...
        output_function_components,...
        output_function_coefficients, glob,params);

descr.output_integral_qdeg = 2; %

%%%%%%%%%% discretization settings

if ~isfield(params,'numintervals_per_block')
  params.numintervals_per_block = 5;
end;

descr.gridtype = 'triagrid';
descr.xnumintervals = params.numintervals_per_block*params.B1;
descr.ynumintervals = params.numintervals_per_block*params.B2;
descr.xrange = [0,1];
descr.yrange = [0,1];

if ~isfield(params,'pdeg')
  params.pdeg = 1; 
end;

if ~isfield(params,'qdeg')
  params.qdeg = 2; 
end;

% numerics settings:
descr.pdeg = params.pdeg; % degree of polynomial functions
descr.qdeg = params.qdeg; % quadrature degree
descr.dimrange = 1; % scalar solution

descr.boundary_type = @thermalblock_boundary_type;
%descr.normals = @my_normals;

% new plot routine, additional block boundaries are plotted
descr.plot_sim_data = @thermalblock_plot_sim_data;
descr.compute_output_functional = 1;
descr.yscale_uicontrols = 0.2;

% for error estimation:
descr.coercivity_alpha = @(descr) min(descr.mus(:));
descr.continuity_gamma = @(descr) max(descr.mus(:));
descr.enable_error_estimator = 1;
descr.get_estimator_from_sim_data = @(sim_data)sim_data.Delta;

% local data functions:
descr = elliptic_discrete_model(descr);

%hp-settings
%------------
descr.distance_function = @euclidean_distance;


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

function model = thermalblock_set_mu(model,mu)
  % updates the parameter
  if length(mu)~=model.number_of_mu
    error('length of mu does not fit to number of blocks!');
  end;
  model.mus = [mu(:)];

function res = thermalblock_boundary_type(glob,params)
  % returns the boundary type
res = zeros(size(glob,1),1);
i  = find(glob(:,1)<=1e-10);
res(i) = -2;
i  = find(glob(:,1)>=1-1e-10);
res(i) = -2;
i  = find(glob(:,2)<=1e-10);
res(i) = -2;
i  = find(glob(:,2)>=1-1e-10);
res(i) = -1;


function comp = thermalblock_neumann_values_components(glob,params) 
  % Neumann values component function
res = zeros(size(glob,1),1);
i = find(glob(:,2)<eps);
res(i) = 1;
comp = {res};

function res = thermalblock_diffusivity_tensor_components(glob,params)
  % diffusivity tensor component function

% 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;
zeroblock = zeros(size(glob,1),4);
res = cell(1,params.number_of_blocks);
for q = 1:params.number_of_blocks;
  block = zeroblock;
  i = find(block_index==q);
  if ~isempty(i)
    block(i,1) = 1;
    block(i,4) = 1;
  end;
  res{q} = block;
end;
% check!
%block = zeroblock;
%for q=1:params.number_of_blocks;
%  block = block + res{q};
%end;
%disp('check diffusivity matrix!')
%keyboard;

function alpha = thermalblock_coercivity_alpha(descr)
  % returns coercivity constant
alpha = min(descr.get_mu(descr));