1 function model = elliptic_debug_model(params);
2 %
function model = elliptic_debug_model(params);
4 % small example of a model, i.e. a structure describing the data
5 % functions and geometry information of a general elliptic equation consisting
6 % of diffusion, convection, reaction equation:
8 % - div ( a(x) grad u(x)) + div (b(x) u(x)) + c(x) u(x) = f(x) on Omega
9 % u(x)) = g_D(x) on Gamma_D
10 % a(x) (grad u(x)) n) = g_N(x) on Gamma_N
11 % alpha(x) u(x) + beta(x) a(x) (grad u(x)) n) = g_R(x) on Gamma_R
13 % s = l(u) linear output functional
15 % Here, we denote the functions as
16 % u: scalar 'solution' (if known, for validation purpose)
18 % a: tensor valued 'diffusivity_tensor'
19 % b: vector valued 'velocity'
20 % c: scalar 'reaction'
22 % g_N: scalar 'neumann_values'
23 % g_R: scalar 'robin_values'
24 % alpha: scalar 'robin_alpha'
25 % beta: scalar 'robin_beta'
27 % Each function allows the evaluation in many points
31 % or model.source(glob,params)
33 % where glob is a n times 2 matrix of row-wise points. The result
34 % is a n times 1 vector of resulting values of f.
36 % model.diffusivity_tensor(glob)
37 % or model.diffusivity_tensor(glob,params)
39 % where glob is a n times 2 matrix of row-wise points. The result
40 % is a n times 4 matrix of resulting values of diffusivity, where
41 % the values of a are sorted in matlab-order as [a_11, a_21, a_12, a_22];
43 % additionally, the model has a function, which determines, whether
44 % a point lies on a Dirichlet, Neumann or Robin boundary:
46 % model.boundary_type(glob)
47 % 0 no boundary (inner edge or point)
48 % -1 indicates Dirichlet-boundary
49 % -2 indicates Neumann-boundary
50 % -3 indicates Robin-boundary
52 % Additionally, the normals in a boundary point can be requested by
55 % Here, glob are assumed to be boundary points lying ON THE
56 % INTERIOR of an edge, such that the outer unit normal is well-defined.
58 % The latter 2 methods boundary_type() and normals() need not be
59 % implemented for grid-based methods, as the normals simply can be
60 % obtained by the grid. The methods are only required, if using
61 % meshless methods with data functions that use normals.
63 % The output functional must be specified
64 % s = model.output_functional(model, model_data, u)
65 % where u is a dof vector of the discrete solution. model_data is
66 % assumed to contain the grid
68 % possible fields of params:
69 % numintervals: the unit square is divided into numintervals
70 % intervals per side. default is 10;
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 % The data functions given in this model are a parametrized model
75 % to be used for debugging purpose. The solution for all parameters
76 % is identical, but the data functions change.
77 % The model parameters are mu_a in [-1,1], mu_b in [0,1] and mu_c
78 % in [0,1]. This model should be used for validating any new
79 % elliptic solver BY PARAMETER VARIATION.
80 % by mu_b=mu_c = 0 the advection and reaction terms can be turned
81 % off. by params.all_dirichlet_boundary = 1, all boundary types can
82 % be set to dirichlet to turn off Neumann and Robin boundaries.
87 % a = [1,0;0,1] + mu_a [0, 1/2; 1/2, 0];
88 % b = [1/2;1/2] (x-y) mu_b;
91 % the exact solution is assumed to be u(x) = sin(k_x x) sin(k_y y)
93 % The computational domain is the unit interval. Dirichlet
94 % boundaries assigned right and upper. Neumann assigned left and
95 % Robin assigned at bottom.
97 % The source term, the boundary values are then specified by the PDE.
99 % An output functional is simply the average over the complete domain
101 % B. Haasdonk 24.1.2011
107 if ~isfield(params,'numintervals')
108 params.numintervals = 10; % 2 x 2 grid! with 8 triangles
111 if ~isfield(params,'pdeg')
115 if ~isfield(params,'qdeg')
120 model = lin_stat_model_default;
124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
125 % Play with the following parameters for debugging a discretization
126 %%% initial setting: simple poisson problem:
130 %model.all_dirichlet_boundary = 1;
131 %%% more complex volume integrals
135 %model.all_dirichlet_boundary = 1;
136 %%% robin and neumann boundary conditions and complex volume integrals
140 model.all_dirichlet_boundary = 0;
142 % parameter settings for rb-methods:
143 model.mu_names ={
'mu_a',
'mu_b',
'mu_c'};
144 model.mu_ranges ={[0,1],[0,1],[0,1]};
145 model.RB_mu_list = {[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],...
146 [1,0,1],[1,1,0],[1,1,1]};
147 model.RB_generation_mode =
'lagrangian';
149 %%%%%%% data
function specification:
151 %
if any of the following flags is set, the corresponding
function
152 % pointers must be provided below.
154 model.has_reaction = 1;
155 %model.has_reaction = 0;
156 model.has_source = 1;
157 model.has_reaction = 1;
158 %model.has_reaction = 0;
159 model.has_advection = 1;
160 %model.has_advection = 0;
161 model.has_diffusivity = 1;
162 model.has_output_functional = 1;
163 model.has_dirichlet_values = 1;
164 model.has_neumann_values = 1;
165 %model.has_neumann_values = 0;
166 %model.has_robin_values = 1;
167 model.has_robin_values = 1;
169 % solution is known:
for validation purpose
170 model.solution = @(glob,params) ...
171 sin(glob(:,1)*params.kx).* ...
172 sin(glob(:,2)*params.ky);
173 % gradient: each row one gradient
174 model.solution_gradient = @(glob,params) ...
175 [params.kx*cos(glob(:,1)*params.kx).* sin(glob(:,2)*params.ky), ...
176 params.ky*sin(glob(:,1)*params.kx).* cos(glob(:,2)*params.ky) ...
178 % hessian: each row: u_xx u_yx, u_xy, u_yy
179 model.solution_hessian = @(glob,params) ...
180 [-params.kx*params.kx*sin(glob(:,1)*params.kx).* sin(glob(:,2)*params.ky), ...
181 params.kx * params.ky * cos(glob(:,1)*params.kx).* cos(glob(:,2)*params.ky) ...
182 params.kx* params.ky * cos(glob(:,1)*params.kx).* cos(glob(:,2)*params.ky), ...
183 -params.ky*params.ky * sin(glob(:,1)*params.kx).* sin(glob(:,2)*params.ky) ...
186 % params is assumed to be the model
187 % c = (1-x) mu_c + x (1-mu_c)
188 reaction_coefficients = @(dummy,params) ...
189 [params.mu_c; 1-params.mu_c];
190 reaction_components = @(glob,params) ...
191 {1-glob(:,1), glob(:,1)};
192 model.reaction = @(glob,params) ...
193 eval_affine_decomp_general(reaction_components, ...
194 reaction_coefficients,glob,params);
196 % b = [1/2; 1/2] *x * mu_b, single component
197 velocity_coefficients = @(dummy,params) ...
199 velocity_components = @(glob,params) ...
200 {ones(size(glob))*0.5.*repmat([glob(:,1)-glob(:,2)],1,2)};
201 model.velocity = @(glob,params) ...
202 eval_affine_decomp_general(velocity_components, ...
203 velocity_coefficients, glob, params);
205 % a = eye(2) + mu_a [0,1/2; 1/2, 0], two components
206 diffusivity_tensor_coefficients = @(dummy,params) ...
208 diffusivity_tensor_components = @(glob,params) ...
210 [ones(size(glob,1),1),...
211 zeros(size(glob,1),1),...
212 zeros(size(glob,1),1),...
213 ones(size(glob,1),1)],...
214 0.5 * [zeros(size(glob,1),1),...
215 ones(size(glob,1),1),...
216 ones(size(glob,1),1),...
217 zeros(size(glob,1),1)] ...
219 model.diffusivity_tensor = @(glob,params) ...
220 eval_affine_decomp_general(diffusivity_tensor_components, ...
221 diffusivity_tensor_coefficients, glob,params);
224 % simply exact solution, not parameter dependent, 1 component!
225 dirichlet_values_coefficients = @(dummy,params) [1];
226 dirichlet_values_components = @(glob,params) ...
227 {model.solution(glob,params)};
228 model.dirichlet_values = @(glob,params) ...
229 eval_affine_decomp_general(dirichlet_values_components, ...
230 dirichlet_values_coefficients,glob,params);
233 % simply exact solution, multiplied with diffusion-tensor and
234 % normal, hence coefficients and component number coincide
235 neumann_values_coefficients = @(dummy,params) ...
236 diffusivity_tensor_coefficients(dummy,params);
237 model.neumann_values = @(glob,params) ...
238 eval_affine_decomp_general(@my_neumann_values_components, ...
239 neumann_values_coefficients, glob, params);
242 % simply weighted sum of dirichlet and neumann values
243 model.robin_alpha = @(glob,params) ones(size(glob,1),1);
244 model.robin_beta = @(glob,params) ones(size(glob,1),1);
245 % first dirichlet coefficients, then neumann coefficients
246 model.robin_values = @(glob,params) ...
247 eval_affine_decomp_general(@my_robin_values_components, ...
248 @my_robin_values_coefficients,glob,params);
250 % source terms: obtain from data functions
251 model.source = @(glob,params) ...
252 eval_affine_decomp_general(@my_source_components, ...
253 @my_source_coefficients,glob,params);
255 % output functional, e.g. average over unit-square:
256 model.output_functional = @output_functional_volume_integral;
257 model.output_function = @output_function_box_mean;
262 model.output_integral_qdeg = 2; %
264 %%%%%%% geometry specification and discretization:
266 model.gridtype =
'triagrid';
267 model.xnumintervals = params.numintervals;
268 model.ynumintervals = params.numintervals;
269 model.xrange = [0,1];
270 model.yrange = [0,1];
273 model.pdeg = params.pdeg; % degree of polynomial functions
274 model.qdeg = params.qdeg; % quadrature degree
275 model.dimrange = 1; % scalar solution
277 % The following 2 methods boundary_type() and normals() need not be
278 % implemented for grid-based methods, as the normals simply can be
279 % obtained by the grid. The methods are only required, if using
280 % meshless methods with data functions that use normals.
281 if model.all_dirichlet_boundary == 1
282 model.boundary_type = @all_dirichlet_boundary_type;
284 model.boundary_type = @mixed_boundary_type;
286 model.normals = @my_normals;
288 % additional settings for reduced basis approach:
289 model.mu_names = {
'mu_a',
'mu_b',
'mu_c'};
290 model.mu_ranges = {[-1,1],[0,1],[0,1]};
292 model.coercivity_alpha = @(model) 1;
293 model.enable_error_estimator = 0;
295 %%%%%%% auxiliary functions:
297 % all edges of unit square are dirichlet, other inner
298 function res = all_dirichlet_boundary_type(glob,params)
299 res = zeros(size(glob,1),1);
300 i = find(glob(:,1)<=1e-10);
301 i = [i, find(glob(:,1)>=1-1e-10)];
302 i = [i, find(glob(:,2)<=1e-10)];
303 i = [i, find(glob(:,2)>=1-1e-10)];
306 % right and upper edges of unit square are dirichlet,
307 % left neumann, lower robin
308 function res = mixed_boundary_type(glob,params)
309 res = zeros(size(glob,1),1);
310 i = find(glob(:,1)<=1e-10);
312 i = find(glob(:,2)<=1e-10);
314 i = find(glob(:,1)>=1-1e-10);
316 i = find(glob(:,2)>=1-1e-10);
319 function res = my_normals(glob,params);
320 res = zeros(size(glob,1),2); % each row one normal
321 i = find(glob(:,1)>1-1e-10);
323 i = find(glob(:,1)<1e-10);
325 i = find(glob(:,2)>1-1e-10);
327 i = find(glob(:,2)<1e-10);
329 % remove diagonal normals
330 i = find (sum(abs(res),2)>1.5);
333 % neumann boundary value components:
334 % multiplication of diffusivity components with solution gradient
335 function res = my_neumann_values_components(glob,params);
336 diff_comp = params.diffusivity_tensor(glob,params);
337 u_grad = params.solution_gradient(glob,params);
338 normals = params.normals(glob,params);
339 res = cell(1,length(diff_comp));
340 for q = 1:length(diff_comp)
341 % Agradu: each row one Agradu value
342 if ~isempty(diff_comp{q})
343 Agradu = [diff_comp{q}(:,1).* u_grad(:,1) + ...
344 diff_comp{q}(:,3).* u_grad(:,2), ...
345 diff_comp{q}(:,2).* u_grad(:,1) + ...
346 diff_comp{q}(:,4).* u_grad(:,2)];
347 res{q} = Agradu(:,1).* normals(:,1) + Agradu(:,2).* normals(:,2);
349 res{q} = zeros(size(glob,1),1);
353 function res = my_robin_values_coefficients(dummy,params)
355 params.neumann_values(dummy,params)];
357 function res = my_robin_values_components(glob,params)
359 neu_comp = params.neumann_values(glob,params);
360 robin_alpha = params.robin_alpha(glob,params);
361 robin_beta = params.robin_beta(glob,params);
362 for q = 1:length(dir_comp)
363 dir_comp{q} = dir_comp{q}.*robin_alpha;
365 for q = 1:length(neu_comp)
366 neu_comp{q} = neu_comp{q}.*robin_beta;
368 res = [dir_comp, neu_comp];
370 % source term: diffusivity, then advection, then reaction contributions
371 function res = my_source_coefficients(dummy,params);
372 res = [-params.diffusivity_tensor(dummy,params);...
373 params.velocity(dummy,params);...
374 params.reaction(dummy,params)];
376 % source term: diffusivity, then advection, then reaction
377 % contribution. Assume that a(x) has vanishing derivative,
378 % i.e piecewise constant!!!
379 % velocity assumed to be divergence free!
380 function res = my_source_components(glob,params);
381 u = params.solution(glob,params);
382 u_grad = params.solution_gradient(glob,params);
383 % hessian: each row: u_xx u_yx, u_xy, u_yy
384 u_hessian = params.solution_hessian(glob,params);
385 diff_comp = params.diffusivity_tensor(glob,params);
386 vel_comp = params.velocity(glob,params);
387 reac_comp = params.reaction(glob,params);
388 source_diff_comp = cell(1,length(diff_comp));
389 source_adv_comp = cell(1,length(vel_comp));
390 source_reac_comp = cell(1,length(reac_comp));
391 for q = 1:length(reac_comp)
392 source_reac_comp{q} = reac_comp{q}.* u;
394 for q = 1:length(vel_comp)
395 source_adv_comp{q} = vel_comp{q}(:,1).* u_grad(:,1) + ...
396 vel_comp{q}(:,2).* u_grad(:,2);
398 for q = 1:length(diff_comp)
399 source_diff_comp{q} = ...
400 diff_comp{q}(:,1).* u_hessian(:,1) + ...
401 diff_comp{q}(:,2).* u_hessian(:,2) + ...
402 diff_comp{q}(:,3).* u_hessian(:,3) + ...
403 diff_comp{q}(:,4).* u_hessian(:,4);
405 res = [source_diff_comp, source_adv_comp, source_reac_comp];
function Udirichlet = dirichlet_values(model, X, Y)
UDIRICHLET = DIRICHLET_VALUES([X],[Y], MODEL) Examples dirichlet_values([0,1,2],[1,1,1],struct(name_dirichlet_values, homogeneous, ... c_dir, 1)) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, xstripes, ... c_dir, [0 1 2], ... dir_borders, [0.3 0.6])) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, box, ... c_dir, 1, ... dir_box_xrange, [0.3 0.6], ... dir_box_yrange, [-0.1 0.1]))