rbmatlab  1.13.10
 All Classes Namespaces Files Functions Variables Groups Pages
demo_femdiscfunc.m
1 function demo_femdiscfunc
2 % function demo_femdiscfunc
3 %
4 % Script demonstrating some basic functionality of lagrange finite
5 % element functions.
6 
7 % B. Haasdonk, I. Maier 26.04.2011
8 
9 disp('---------------------------------');
10 disp(' lagrange FE-functions ');
11 disp('---------------------------------');
12 
13 % poisson_model
14 params = [];
15 pdeg = 4;
16 params.pdeg = pdeg;
17 params.dimrange = 3;
18 % params.dimrange = 1;
19 params.debug = 1;
20 params.numintervals = 5;
21 model = poisson_model(params);
22 % convert to local_model:
23 model = elliptic_discrete_model(model);
24 grid = construct_grid(model);
25 df_info = feminfo(params,grid);
26 %tmp = load('circle_grid');
27 %grid = triagrid(tmp.p,tmp.t,[]);
28 disp('model initialized');
29 
30 % without model_data, detailed_simulation, etc. but explicit
31 % calls of fem operations
32 
33 % initialize vectorial discrete function, extract scalar
34 % component and plot basis function
35 df = femdiscfunc([],df_info); % initialize zero df
36 
37 fprintf('\n');
38 disp('initialization of femdiscfunc successful, display:');
39 display(df);
40 
41 disp('press enter to continue');
42 pause;
43 
44 % for check of dof consistency we have a plot routine:
45 fprintf('\n');
46 disp('plot of global_dof_index map for consistency check:');
47 p = plot_dofmap(df);
48 
49 disp('press enter to continue');
50 pause;
51 
52 % global dof access:
53 fprintf('\n');
54 disp('example of dof access, scalar component extraction and plot:');
55 % set second vector component in 4th basis function nonzero;
56 ncomp = 2;
57 locbasisfunc_index = 4;
58 df.dofs((locbasisfunc_index-1)*df.dimrange+ncomp) = 1;
59 dfscalar = scalar_component(df,2); % should be nonzero function
60 disp(['entry should be 1 : ',...
61  num2str(dfscalar.dofs(locbasisfunc_index))]);
62 
63 params = [];
64 params.subsampling_level = 6;
65 figure,plot(dfscalar,params);
66 dfscalar.dofs(:) = 0; % now zero again.
67 % keyboard;
68 
69 disp('press enter to continue');
70 pause;
71 
72 fprintf('\n');
73 disp('example of local dof access:');
74 % local dof access via local to global map
75 eind = 4; % set dof on element number 4
76 lind = (pdeg+2); % local basis function index with lcoord (0,1/pdeg)
77 dfscalar.dofs(dfscalar.global_dof_index(eind,lind)) = 1;
78 
79 % example of local evaluation of femdiscfunc simultaneous on
80 % several elements in the same local coordinate point
81 elids = [4,6,7,10]; % evaluation on some elements
82 lcoord = [0,1/pdeg]; % local coordinate vector == last point in all triangles
83 f = evaluate(dfscalar,elids,lcoord);
84 disp(['first entry should be 1 : ',num2str(f(1))]);
85 % equivalent call (!) by () operator as abbreviation for local evaluation:
86 f = dfscalar(elids,lcoord);
87 disp(['first entry should be 1 : ',num2str(f(1))]);
88 
89 disp('press enter to continue');
90 pause;
91 
92 disp('examples of norm computation:')
93 params.dimrange = 1;
94 params.pdeg = 1;
95 dfinfo1 = feminfo(params,grid);
96 df1 = femdiscfunc([],dfinfo1);
97 df1.dofs(:) = 1;
98 disp(['L2-norm(f(x,y)=1) = ',num2str(fem_l2_norm(df1))]);
99 disp(['H10-norm(f(x,y)=1) = ',num2str(fem_h10_norm(df1))]);
100 df1.dofs(:) = df1.grid.X(:);
101 disp(['L2-norm(f(x,y)=x) = ',num2str(fem_l2_norm(df1))]);
102 disp(['H10-norm(f(x,y)=x) = ',num2str(fem_h10_norm(df1))]);
103 df1.dofs(:) = df1.grid.Y(:);
104 disp(['L2-norm(f(x,y)=y) = ',num2str(fem_l2_norm(df1))]);
105 disp(['H10-norm(f(x,y)=y) = ',num2str(fem_h10_norm(df1))]);
106 
107 disp('press enter to continue');
108 pause;
109 
110 % evaluate df in all lagrange nodes of element 4 by loop
111 fprintf('\n');
112 disp(['dfscalar on element 4 in all lagrange nodes,' ...
113  'only (pdeg+2) entry should be 1:']);
114 lagrange_nodes = lagrange_nodes_lcoord(pdeg);
115 elid = 4;
116 for i = 1:size(lagrange_nodes,1);
117  f = evaluate(dfscalar,elid,lagrange_nodes(i,:));
118  disp(['f(l(',num2str(i),')) = ',num2str(f)]);
119 end;
120 
121 disp('press enter to continue');
122 pause;
123 
124 fprintf('\n');
125 disp('example of requirement of subsampling in plot of discfuncs:');
126 figure;
127 subsamp_levels = [0,2,4,16];
128 for i=1:length(subsamp_levels)
129  subplot(2,2,i),
130  params.axis_equal = 1;
131  params.subsampling_level = subsamp_levels(i);
132  params.clim = [-0.15 1.15]; % rough bounds
133  plot(dfscalar,params);
134  title(['subsampling level = ',num2str(subsamp_levels(i))]);
135 end;