1 function model = poisson_model(params);
2 %
function model = poisson_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'
21 % g_D: scalar 'dirichlet_values'
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 the simple poisson
75 % equation with Gamma_D = boundary(Omega), Gamma_N = {}, Gamma_R = {}
80 % with exact solution u(x) = 16 x_1(1-x_1)x_2(1-x_2),
81 % i.e. f(x) = 32 (x_1 + x_2 - x_1^2 - x_2^2)
83 % B. Haasdonk 23.11.2010
89 if ~isfield(params,
'numintervals')
90 params.numintervals = 10; % 2 x 2 grid! with 8 triangles
93 if ~isfield(params,'pdeg')
94 params.pdeg = 2; % 2 x 2 grid! with 8 triangles
97 if ~isfield(params,'qdeg')
98 params.qdeg = 2; % 2 x 2 grid! with 8 triangles
102 model.rb_problem_type = 'lin_elliptic';
104 %%%%%%% data function specification:
106 % if any of the following flags is set, the corresponding function
107 % pointers must be provided below.
109 model.has_reaction = 1;
110 %model.has_reaction = 0;
111 model.has_source = 1;
112 model.has_reaction = 1;
113 %model.has_reaction = 0;
114 model.has_advection = 1;
115 %model.has_advection = 0;
116 model.has_diffusivity = 1;
117 model.has_output_functional = 1;
118 model.has_dirichlet_values = 1;
119 model.has_neumann_values = 1;
120 %model.has_neumann_values = 0;
121 %model.has_robin_values = 1;
122 model.has_robin_values = 1;
125 model.solution = @(glob,params) ...
126 16 * glob(:,1).*(1-glob(:,1)).*glob(:,2).*(1-glob(:,2));
127 % for debugging: f = 1; u = -1/4 (x^2+y^2)
128 %model.solution = @(glob,params) ...
129 % -0.25 * (glob(:,1).^2+glob(:,2).^2);
131 % params is an optional parameter, perhaps useful later
132 model.source = @(glob,params) ...
133 32 * (glob(:,1)+ glob(:,2)- glob(:,1).^2 - glob(:,2).^2);
134 % for debugging: f = 1; u = -0.25* (x^2+y^2)
135 %model.source = @(glob,params) ...
136 % ones(size(glob,1),1);
138 model.reaction = @(glob,params) zeros(size(glob,1),1);
139 model.velocity = @(glob,params) zeros(size(glob,1),2);
140 model.diffusivity_tensor = @(glob,params) ...
141 [ones(size(glob,1),1),...
142 zeros(size(glob,1),1),...
143 zeros(size(glob,1),1),...
144 ones(size(glob,1),1)];
145 %model.diffusivity_gradient = @(glob,params) zeros(size(glob,1),2);
146 model.reaction = @(glob,params) zeros(size(glob,1),1);
148 % Dirichlet everywhere, see function below.
151 params.solution(glob,params);
152 model.neumann_values = @(glob,params) zeros(size(glob,1),1);
153 model.robin_values = @(glob,params) zeros(size(glob,1),1);
154 model.robin_alpha = @(glob,params) zeros(size(glob,1),1);
155 model.robin_beta = @(glob,params) zeros(size(glob,1),1);
157 % output functional, e.g. average over unit-square:
158 model.output_functional = @output_functional_volume_integral;
159 model.output_function = @output_function_box_mean;
164 model.output_integral_qdeg = 0; % midpoint integration
165 % later, e.g. for thermal block:
166 %model.output_functional = @output_functional_boundary_integral;
167 %model.output_function = @output_function_dirichlet_indicator;
169 %%%%%%% geometry specification and discretization:
172 model.xnumintervals = params.numintervals;
173 model.ynumintervals = params.numintervals;
174 model.xrange = [0,1];
175 model.yrange = [0,1];
178 model.pdeg = params.pdeg; % degree of polynomial functions
179 model.qdeg = params.qdeg; % quadrature degree
180 model.dimrange = 1; % scalar solution
182 % The following 2 methods boundary_type() and normals() need not be
183 % implemented for grid-based methods, as the normals simply can be
184 % obtained by the grid. The methods are only required, if using
185 % meshless methods with data functions that use normals.
186 model.boundary_type = @my_boundary_type;
187 model.normals = @my_normals;
189 %%%%%%% auxiliary functions:
191 % all edges of unit square are dirichlet, other inner
192 % note this function is not used by grid-based methods
193 function res = my_boundary_type(glob,params)
194 res = zeros(size(glob,1),1);
195 i = find(glob(:,1)<=1e-10);
196 i = [i, find(glob(:,1)>=1-1e-10)];
197 i = [i, find(glob(:,2)<=1e-10)];
198 i = [i, find(glob(:,2)>=1-1e-10)];
201 function res = my_normals(glob,params);
202 res = zeros(size(glob,1),2); % each row one normal
203 i = find(glob(:,1)>1-1e-10);
205 i = find(glob(:,1)<-1+1e-10);
207 i = find(glob(:,2)>1-1e-10);
209 i = find(glob(:,2)<-1+1e-10);
211 % remove diagonal normals
212 i = find (sum(abs(res),2)>1.5);