2 % Thermal Block example similar as described in the book of
3 % A.T. patera and G. Rozza (just one parameter more)
6 % - div ( a(x) grad u(x)) = q(x) on Omega
7 % u(x)) = g_D(x) on Gamma_D
8 % a(x) (grad u(x)) n) = g_N(x) on Gamma_N
9 % s = l(u) linear output functional
11 % where Omega = [0,1]^2 is divided into B1 * B2 blocks
12 % QA := B1*B2. The heat conductivities are given by mu_i:
14 % -------------------------
15 % | ... | ..... | mu_QA |
16 % -------------------------
17 % | .... | ... | ..... |
18 % -------------------------
19 % | mu_1 | ... | mu_B1 |
20 % -------------------------
22 % `Gamma_D =` upper edge
23 % `Gamma_N = boundary(Omega)\Gamma_D`
25 % `a(x) = mu_i(x) if x\in B_i`
28 % `g_N = 1` on lower edge, 0 otherwise
29 % `l(u)` = average over lower edge
31 % S. Langhof and B. Haasdonk 26.2.2011
34 descr = lin_stat_model_default;
35 descr.name =
'thermalblock';
36 descr.dmodel_constructor = @ThermalBlock.DetailedModel;
38 descr.set_mu = @thermalblock_set_mu;
39 descr.get_mu = @(descr) descr.mus(:);
44 % descr.default =
'3x3 demo version';
46 % number of blocks in X and Y direction
49 descr.number_of_blocks = params.B1*params.B2;
51 % each block has a different heat conductivity
52 % upper left has conductivity 1
56 for p = 1:descr.number_of_blocks
57 mu_names = [mu_names,{['mu',num2str(p)]}];
58 mu_ranges = [mu_ranges,{mu_range}];
60 descr.mu_names = mu_names;
61 descr.mu_ranges = mu_ranges;
63 %
default values 1 everywhere
64 descr.mus = ones(descr.number_of_blocks,1);
66 % simple snapshots without orthogonalization
67 descr.RB_mu_list = {[0.1,1,1,1,1,1,1,1,1],...
68 [1,0.1,1,1,1,1,1,1,1],...
69 [1,1,0.1,1,1,1,1,1,1],...
70 [1,1,1,0.1,1,1,1,1,1],...
71 [1,1,1,1,0.1,1,1,1,1],...
72 [1,1,1,1,1,0.1,1,1,1],...
73 [1,1,1,1,1,1,0.1,1,1],...
74 [1,1,1,1,1,1,1,0.1,1],...
75 [1,1,1,1,1,1,1,1,0.1]};
77 %%%%%%%%%
set data functions
78 descr.has_diffusivity = 1;
79 descr.has_output_functional = 1;
80 descr.has_dirichlet_values = 1;
81 descr.has_neumann_values = 1;
83 % zero dirichlet values, i.e. 1 component, Q_dir = 1
84 dirichlet_values_coefficients = @(dummy,params) [0];
85 dirichlet_values_components = @(glob,params) {zeros(size(glob,1),1)};
86 descr.dirichlet_values = @(glob,params) ...
87 eval_affine_decomp_general(dirichlet_values_components, ...
88 dirichlet_values_coefficients,glob,params);
90 % 1/0 neumann values depending on edge, non parametric, i.e Q_neu = 1;
91 neumann_values_coefficients = @(dummy,params) 1; % single component
92 descr.neumann_values = @(glob,params) ...
93 eval_affine_decomp_general(@thermalblock_neumann_values_components, ...
94 neumann_values_coefficients, glob, params);
96 % diffusion tensor: each row four entries a11,a_21,a_12,a_22.
97 % a11(x)=a22(x) = mu_i
if x in block i, a12=a21 = 0.
98 diffusivity_tensor_coefficients = @(dummy,params) ...
100 descr.diffusivity_tensor = @(glob,params) ...
101 eval_affine_decomp_general(...
102 @thermalblock_diffusivity_tensor_components, ...
103 diffusivity_tensor_coefficients, glob,params);
105 % only useful
for detailed simulation or nonlinear outputs:
106 %descr.output_functional = @output_functional_boundary_integral;
108 descr.operators_output = @fem_operators_output_boundary_integral;
109 % output weight
function: simply 1 on lower edge, 0
else => Ql = 1 component
110 output_function_components = @(glob,descr) {1.0*(glob(:,2)<eps)};
111 output_function_coefficients = @(glob,descr) 1;
112 descr.output_function = @(glob,params) ...
113 eval_affine_decomp_general(...
114 output_function_components,...
115 output_function_coefficients, glob,params);
117 descr.output_integral_qdeg = 2; %
119 %%%%%%%%%% discretization settings
121 if ~isfield(params,
'numintervals_per_block')
122 params.numintervals_per_block = 5;
126 descr.xnumintervals = params.numintervals_per_block*params.B1;
127 descr.ynumintervals = params.numintervals_per_block*params.B2;
128 descr.xrange = [0,1];
129 descr.yrange = [0,1];
131 if ~isfield(params,'pdeg')
135 if ~isfield(params,'qdeg')
140 descr.pdeg = params.pdeg; % degree of polynomial functions
141 descr.qdeg = params.qdeg; % quadrature degree
142 descr.dimrange = 1; % scalar solution
144 descr.boundary_type = @thermalblock_boundary_type;
145 %descr.normals = @my_normals;
147 % new plot routine, additional block boundaries are plotted
149 descr.compute_output_functional = 1;
150 descr.yscale_uicontrols = 0.2;
152 % for error estimation:
153 descr.coercivity_alpha = @(descr) min(descr.mus(:));
154 descr.continuity_gamma = @(descr) max(descr.mus(:));
155 descr.enable_error_estimator = 1;
156 descr.get_estimator_from_sim_data = @(sim_data)sim_data.Delta;
158 % local data functions:
159 descr = elliptic_discrete_model(descr);
161 %%%%%%% auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163 function model = thermalblock_set_mu(model,mu)
164 % updates the parameter
165 if length(mu)~=model.number_of_blocks
166 error('length of mu does not fit to number of blocks!');
170 function res = thermalblock_boundary_type(glob,params)
171 % returns the boundary type
172 res = zeros(size(glob,1),1);
173 i = find(glob(:,1)<=1e-10);
175 i = find(glob(:,1)>=1-1e-10);
177 i = find(glob(:,2)<=1e-10);
179 i = find(glob(:,2)>=1-1e-10);
183 function comp = thermalblock_neumann_values_components(glob,params)
184 % Neumann values component function
185 res = zeros(size(glob,1),1);
186 i = find(glob(:,2)<eps);
190 function res = thermalblock_diffusivity_tensor_components(glob,params)
191 % diffusivity tensor component
function
193 %
for each point in glob find global block number
194 % xi: range 0 ... B1-1
195 xi = floor(glob(:,1)*params.B1);
196 i = find(xi>=params.B1);
200 % xi: range 0 ... B1-1
201 yi = floor(glob(:,2)*params.B2);
202 i = find(yi>=params.B2);
206 block_index = yi*params.B1+xi+1;
207 zeroblock = zeros(size(glob,1),4);
208 res = cell(1,params.number_of_blocks);
209 for q = 1:params.number_of_blocks;
211 i = find(block_index==q);
220 %
for q=1:params.number_of_blocks;
221 % block = block + res{q};
223 %disp(
'check diffusivity matrix!')
226 function alpha = thermalblock_coercivity_alpha(descr)
227 % returns coercivity constant
228 alpha = min(descr.get_mu(descr));