1 function res = fem_poisson(step,params)
2 %
function res = fem_poisson(step)
4 % Script demonstrating lagrange finite element functions,
using
5 % them
for function interpolation and finite element methods.
7 % Script realizing the simple poisson equation or a more complex
8 % elliptic problem with arbitrary diffusivity_tensor, advection,
9 % source, reaction, neumann , dirichlet and robin boundary
10 % conditions on the unit square
11 % with finite elements on
triagrid discretization.
13 % step 1: some basic
femdiscfunc demonstration and functionatity
14 % step 2: use
femdiscfunc for interpolating functions by local and
16 % step 3: solve fem-problem, return errors
17 % step 4 error convergence investigation
18 % step 5: same as 3 by call of interface methods + Reduced basis steps
21 % B. Haasdonk 11.1.2011
23 % adapted to new assembly
25 % Immanuel Maier 12.04.2011
41 case 1 % elementary fem discrete
function operations
47 % params.dimrange = 1;
49 params.numintervals = 5;
50 model = poisson_model(params);
51 % convert to local_model:
52 model = elliptic_discrete_model(model);
53 grid = construct_grid(model);
55 %tmp = load(
'circle_grid');
57 disp(
'model initialized');
59 % without model_data, detailed_simulation, etc. but
explicit
60 % calls of fem operations
62 % initialize vectorial discrete
function, extract scalar
63 % component and plot basis
function
67 disp(
'initialization of femdiscfunc successful, display:');
70 disp(
'press enter to continue');
73 %
for check of dof consistency we have a plot routine:
75 disp(
'plot of global_dof_index map for consistency check:');
78 disp(
'press enter to continue');
83 disp(
'example of dof access, scalar component extraction and plot:');
84 %
set second vector component in 4th basis
function nonzero;
86 locbasisfunc_index = 4;
87 df.dofs((locbasisfunc_index-1)*df.dimrange+ncomp) = 1;
88 dfscalar = scalar_component(df,2); % should be nonzero
function
89 disp([
'entry should be 1 : ',...
90 num2str(dfscalar.dofs(locbasisfunc_index))]);
93 params.subsampling_level = 6;
94 figure,plot(dfscalar,params);
95 dfscalar.dofs(:) = 0; % now zero again.
98 disp(
'press enter to continue');
102 disp(
'example of local dof access:');
103 % local dof access via local to global map
104 eind = 4; %
set dof on element number 4
105 lind = (pdeg+2); % local basis
function index with lcoord (0,1/pdeg)
106 dfscalar.dofs(dfscalar.global_dof_index(eind,lind)) = 1;
108 % example of local evaluation of
femdiscfunc simultaneous on
109 % several elements in the same local coordinate point
110 elids = [4,6,7,10]; % evaluation on some elements
111 lcoord = [0,1/pdeg]; % local coordinate vector == last point in all triangles
112 f = evaluate(dfscalar,elids,lcoord);
113 disp(['first entry should be 1 : ',num2str(f(1))]);
114 % equivalent call (!) by () operator as abbreviation for local evaluation:
115 f = dfscalar(elids,lcoord);
116 disp(['first entry should be 1 : ',num2str(f(1))]);
118 disp('press enter to continue');
121 disp('examples of norm computation:')
124 dfinfo1 =
feminfo(params,grid);
127 disp(['L2-norm(f(x,y)=1) = ',num2str(fem_l2_norm(df1))]);
128 disp(['H10-norm(f(x,y)=1) = ',num2str(fem_h10_norm(df1))]);
129 df1.dofs(:) = df1.grid.X(:);
130 disp(['L2-norm(f(x,y)=x) = ',num2str(fem_l2_norm(df1))]);
131 disp(['H10-norm(f(x,y)=x) = ',num2str(fem_h10_norm(df1))]);
132 df1.dofs(:) = df1.grid.Y(:);
133 disp(['L2-norm(f(x,y)=y) = ',num2str(fem_l2_norm(df1))]);
134 disp(['H10-norm(f(x,y)=y) = ',num2str(fem_h10_norm(df1))]);
136 disp('press enter to continue');
139 % evaluate df in all lagrange nodes of element 4 by loop
141 disp(['dfscalar on element 4 in all lagrange nodes,' ...
142 'only (pdeg+2) entry should be 1:']);
143 lagrange_nodes = lagrange_nodes_lcoord(pdeg);
145 for i = 1:size(lagrange_nodes,1);
146 f = evaluate(dfscalar,elid,lagrange_nodes(i,:));
147 disp(['f(l(',num2str(i),')) = ',num2str(f)]);
150 disp('press enter to continue');
154 disp('example of requirement of subsampling in plot of discfuncs:');
156 subsamp_levels = [0,2,4,16];
157 for i=1:length(subsamp_levels)
159 params.axis_equal = 1;
160 params.subsampling_level = subsamp_levels(i);
161 params.clim = [-0.15 1.15]; % rough bounds
162 plot(dfscalar,params);
163 title(['subsampling level = ',num2str(subsamp_levels(i))]);
166 disp('end of step 1, please inspect variables, if required.')
175 params.numintervals = 5;
176 model = poisson_model(params);
177 % convert to local_model:
178 model = elliptic_discrete_model(model);
179 grid = construct_grid(model);
180 %tmp = load('circle_grid');
182 disp('model initialized');
184 % interpolate exact solution and other coefficient functions
185 % as fem-function and plot
187 disp('examples of interpolation of analytical functions:');
189 df = femdiscfunc([],df_info);
191 df = fem_interpol_global(model.solution, df);
193 title('analytical solution');
195 % discretize source function and plot
197 % problems with fem_interpol_local!!
198 %df = fem_interpol_local(model.source, df);
199 df = fem_interpol_global(model.source, df);
201 title('source function');
203 disp('press enter to continue');
206 % discretize diffusivity and plot as 4-sequence
207 % to be implemented for vectorial functions...
208 %params.dimrange = 4;
209 %df4 = femdiscfunc([],grid,params);
210 %df4 = fem_interpol_global(model.diffusivity_tensor, df)
211 % arrange as sequence of 4 scalar functions and plot
213 %title = 'diffusivity function';
215 % example of local evaluation of vectorial basis function
216 % on reference triangle
217 disp('local evaluation of vectorial basis functions on reference triangle:');
222 df2 = femdiscfunc([],df_info);
223 res = fem_evaluate_basis(df2,lcoord)
225 disp('local evaluation of scalar basis functions derivative on reference triangle:');
226 res = fem_evaluate_scalar_basis_derivative(df,lcoord)
229 % disp('local evaluation of vectorial basis function derivative on reference triangle:');
232 % interpolation error convergence with grid refinement
234 case 3 % finite element assembly, solution and plot.
235 % different components assembled separately, as later
236 % affine-decomposition will require such modularization
238 % assemble discrete system, solve and plot results
244 qdeg = 3 * pdeg; % quadrature_degree
247 params.numintervals = 20;
248 params.no_plot = 0; %plotting active
251 no_plot = params.no_plot;
253 %model = poisson_model(params);
254 model = elliptic_debug_model(params);
255 model = elliptic_discrete_model(model); % convert to local_model
256 grid = construct_grid(model);
257 disp(['grid.nelements=',num2str(grid.nelements)]);
258 % for debugging: set all to neumann
259 % i = find(grid.NBI==-1);
263 disp(['df_info.ndofs=',num2str(df_info.ndofs)]);
264 %tmp = load('circle_grid');
266 disp('model initialized');
268 %%%%%%%%%%%%%%% assembly of system %%%%%%%%%%%%%%%%%%%
270 % assemble right hand side
271 [r_source, r_dirichlet, r_neumann, r_robin] = ...
272 fem_rhs_parts_assembly(model,df_info);
273 r = r_source + r_neumann + r_robin + r_dirichlet;
274 disp('rhs assembled');
276 % sparse system matrix:
277 A = spalloc(df_info.ndofs,df_info.ndofs,10); % rough upper bound for number of nonzeros
278 [A_diff , A_adv, A_reac, A_dirichlet, A_neumann, A_robin] = ...
279 fem_matrix_parts_assembly(model,df_info);
280 A = A_diff + A_adv + A_reac + A_neumann + A_robin + A_dirichlet;
281 disp('matrix assembled');
283 %%%%%%%%%%%%%%% solve system %%%%%%%%%%%%%%%%%%%
285 u_h = femdiscfunc([],df_info);
287 disp('system solved');
289 %%%%%%%%%%%%%%% postprocessing %%%%%%%%%%%%%%%%%%%
294 params.subsampling_level = 10;
295 figure, plot(u_h,params);
296 title('discrete solution')
299 df = femdiscfunc([],df_info);
300 df = fem_interpol_global(model.solution, df,model);
303 title('analytical solution');
306 % plot error of (interpolated) exact solution and discrete sol.
307 % interpolation using same degree as discrete solution
308 e_h = femdiscfunc([],df_info);
309 e_h = fem_interpol_global(model.solution, e_h,model);
310 e_h.dofs = e_h.dofs -u_h.dofs;
313 title('rough error I\_rough(u)-u_h');
317 % better: choose high degree for interpolation of u
318 % then require local_interpolation of u_h into high-resolved space
319 params.pdeg = 4; % choose high degree
321 % map u_h to higher degree discretefunction
322 params.has_dirichlet_values = 10;
323 df_info_fine =
feminfo(params,grid);
324 u_h_fine = femdiscfunc([],df_info_fine);
325 u_h_fine = fem_interpol_local(u_h, u_h_fine,model);
326 e_h_fine = femdiscfunc([],df_info_fine);
327 e_h_fine = fem_interpol_global(model.solution, e_h_fine,model);
328 e_h_fine.dofs = e_h_fine.dofs -u_h_fine.dofs;
329 figure, plot(e_h_fine);
330 title('fine error I\_fine(u)-u_h');
333 % determine error in L2 and H10 norm
334 l2_err = fem_l2_norm(e_h);
335 h10_err = fem_h10_norm(e_h);
336 disp(['L2-error norm |u_h - I(u_exact)| = ',num2str(l2_err)]);
337 disp(['H10-error norm |u_h - I(u_exact)| = ',num2str(h10_err)]);
340 disp('solution and errors plotted. Inspect workspace.')
344 res.l2_error = l2_err;
345 res.h10_error = h10_err;
347 case 4 % error convergence investigation
352 qdeg = 3 * pdeg; % quadrature_degree
357 numintervals = [2,4,8,16,32,64,128,256];
358 l2_errs = zeros(length(numintervals),1);
359 h10_errs = zeros(length(numintervals),1);
360 for i = 1:length(numintervals)
361 disp(['----------------------------------------------']);
362 disp(['numintervals = ',num2str(numintervals(i))]);
363 params.numintervals = numintervals(i);
364 res = fem_poisson(3,params);
365 l2_errs(i) = res.l2_error;
366 h10_errs(i) = res.h10_error;
369 case 5 % detailed and rb simulation by high-level methods
371 disp('detailed simulation:');
372 model = elliptic_debug_model;
373 model = elliptic_discrete_model(model);
374 model_data = gen_model_data(model);
375 model = model.set_mu(model,[0,0,0]);
376 sim_data = detailed_simulation(model,model_data);
377 plot_params.title = 'detailed_solution';
378 figure,
plot_sim_data(model,model_data,sim_data,plot_params);
380 disp('offline computations (gen_detailed_data: basis):');
381 detailed_data = gen_detailed_data(model,model_data);
382 disp('offline computations II (gen_reduced_data: operator components):');
383 reduced_data = gen_reduced_data(model,detailed_data);
385 disp('online reduced simulation:');
386 model = model.set_mu(model,[0,0,0]);
387 rb_sim_data = rb_simulation(model,reduced_data);
388 rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
389 plot_params.title = 'reduced solution';
390 figure,
plot_sim_data(model,model_data,rb_sim_data,plot_params);
392 eh = sim_data.uh - rb_sim_data.uh;
393 l2_err = fem_l2_norm(eh);
394 h10_err = fem_h10_norm(eh);
395 disp(['L2-error: ',num2str(l2_err)]);
396 disp(['H10-error: ',num2str(h10_err)]);
401 model = elliptic_debug_model;
402 model = elliptic_discrete_model(model);
403 model_data = gen_model_data(model);
404 detailed_data = gen_detailed_data(model,model_data);
405 plot_params.axis_tight = 1;