3 % demo
for showing ldgfunc capabilities
5 % idea is mainly to keep dofs and size information separate in
6 % order to easily work with vectors and matrices of dofs
for
9 % A slower alternative can be beneficial: by
using the \\@ldg
function
10 %
class, a dof-vector can be stored in the
object. And hence an
11 % evaluation be performed by the subsref method as
for analytical functions.
12 % So ldg-objects can be used identically as any
13 % analytical
function in integration, matrix assembly, etc.
14 % As classes/methods are slower than
15 % structs,
this should only be used
if necessary.
16 % So most ldg operations will accept separate dof vectors/matrices
17 % and size information.
19 % Bernard Haasdonk 28.1.2009
21 % quadrature degree
for integration:
27 disp(
'initialization of zero ldg function');
29 % initialize basic data
36 df_info =
ldginfo(params, grid);
37 %params.nelements = grid.nelements;
38 %params.ndofs = ldg_ndofs(params);
39 %params.ndofs_per_element = ldg_ndofs_per_element(params);
40 dofs = ldg_zero(df_info); % simple zero vector of dofs
42 % these variables are sufficient to perform most effective operations
43 %
for ldg functions, i.e. data and parameters are kept disjoint.
45 % A slower alternative can be beneficial: by
using the @ldg
function
46 %
class, a dof-vector can be stored in the
object. And hence an
47 % evaluation be performed by the subsref method as
for analytical functions.
48 % So ldg-objects can be used identically as any
49 % analytical
function in integration, matrix assembly, etc.
50 % As classes/methods are slower than
51 % structs,
this should only be used
if necessary:
53 df =
ldgdiscfunc(dofs,df_info); % setting of data: persistent!
56 % most of the following does not access df, but work with dofs and params
58 disp(
'----------------------------------------------------------');
59 disp([
'dof-access and plot of scalar component with subsampling']);
61 % manipulate/read dofs and visualize
62 df.dofs = 1:df_info.ndofs;
63 [dofs_scalar1,sparams] = ldg_scalar_component(df,1);
67 plot_params.plot_subsampling_level = ssl;
69 ldg_plot(df_scalar,grid,plot_params);
70 title([
'dofnumbers, subsampling level ',num2str(ssl)]);
72 disp(
'press key to continue');
76 disp(
'----------------------------------------------------------');
77 disp([
'local evaluation of ldgfunc and analytical functions in point ',...
78 ' lcoord=[1/3,1/3] in first 10 elements:']);
79 % define analytical
function inline, e.g. constant 1
80 f1 = @(einds,loc,grid,params) ones(length(einds),1);
81 % define analytical
function by pointer to matlab
function
82 f2 = @(einds,loc,grid,params) f_local(einds,loc,grid,params);
83 % define analytical
function by pointer to matlab
function with
84 %
internal conversion of local to global coordinates:
85 f3 = @(einds,loc,grid,params) f_global(...
86 local2global(grid,einds,loc,params),params);
87 % evaluate all functions
88 disp(
'f1 (scalar function constant 1):');
89 f1(1:10, [1/3,1/3], grid, params)
90 disp('f2: (function f_local defined in local coordinates)');
91 f2(1:10, [1/3,1/3], grid, params)
92 disp('f3: (function f_global defined in global coordinates)');
93 f3(1:10, [1/3,1/3], grid, params)
94 % note here a ldg function is used identically as the analytical functions!
95 disp('df: (@
ldgdiscfunc dimrange=2, pdeg=1, zero-function)');
96 df.dofs = ldg_zero(df_info);
98 df(1:10, [1/3,1/3]) % identical call as above!
99 % this is equivalent to
100 % ldg_evaluate(dofs,1:10, [1/3,1/3], grid, params);
102 disp('press key to continue');
106 disp('----------------------------------------------------------');
107 disp('test of orthogonality of element basis functions: gram matrix');
109 % test orthogonality of basis functions:
110 disp('gram matrix of basis on reference element should be identity:')
111 K = ldg_local_mass_matrix(qdeg,params)
112 %f = @(lcoord) gram_matrix(ldg_evaluate_basis(lcoord,params));
113 %triaquadrature(qdeg,f)
115 disp('press key to continue');
119 disp('----------------------------------------------------------');
120 disp('l2 projection of analytic function: hor/vert sinus waves');
122 % l2 projection of an analytical function and plot
123 f = @(einds,loc,grid,varargin) f_global(...
124 local2global(grid,einds,loc,df_info),df_info);
125 df.dofs = ldg_l2project(f,qdeg,df_info);
126 %dofs = l2project(f,qdeg,grid,params);
128 [dofs_scalar1,sparams] = ldg_scalar_component(df,1);
129 dofs_scalar2 = ldg_scalar_component(df,2);
132 plot_params.plot_subsampling_level = 2;
133 figure,ldg_plot(df_scalar1,grid,plot_params);title('component1');
134 figure,ldg_plot(df_scalar2,grid,plot_params);title('component2');
136 disp('press key to continue');
140 disp('----------------------------------------------------------');
141 disp('l2-error of zero ldg function with unity (should be 1)');
146 params = ldg_params(grid.nelements,pdeg,dimrange);
147 df_info =
ldginfo(params, grid);
149 % initialize zero function:
151 df.dofs = ldg_zero(params);
153 % initialize const-one function analytically
154 f1 = @(einds,loc,grid,params) ones(length(einds),1);
155 % and project to ldgfunc
156 dofs1 = ldg_l2project(f1,qdeg,df_info);
158 df1 = ldgdiscfunc(dofs1,df_info);
161 % compute errors with analytical and discrete function:
162 % identical syntax by using ldg class!!!
163 params.evaluate = @ldg_evaluate;
167 err_df0_f1 = ldg_l2error(df0,f1,qdeg);
168 disp(['l2 error (df0,f1) =',num2str(err_df0_f1)]);
169 err_df0_df1 = ldg_l2error(df0,df1,qdeg);
170 disp(['l2 error (df0,df1) =',num2str(err_df0_df1)]);
172 disp('press key to continue');
176 disp('----------------------------------------------------------');
177 disp(['Table of l2-errors of projection with increasing p, should' ...
179 % l2-difference between ldgdiscfunc and function with
180 % entity-local-coord evaluation possibility.
184 params = ldg_params(grid.nelements,p,dimrange);
185 df_info =
ldginfo(params, grid);
187 dofs = ldg_l2project(f,qdeg,df_info);
188 df = ldgdiscfunc(dofs,df_info);
190 l2err = ldg_l2error(df,f,qdeg);
191 disp(['pol-deg = ',num2str(p),', l2error(df,f) = ',num2str(l2err)]);
194 disp('press key to continue');
199 % l2-difference between two ldgdiscfuncs works identical
200 disp(['l2error(df,df) =',...
201 num2str(l2error(df,df,qdeg,grid,params))]);
203 % arithmetics of ldg functions: sum, diff, scal-mtimes,
204 % in particular scalar product of vectorial localdgs
208 % right hand side assembly of LGS
220 % div0 projection operator
228 % concatenate for complete space operator
233 % implicit or explicit solution of time step
238 % stiffness-matrix assembly into sparse matrix
245 % stop routine for inspection of workspace variables
248 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
249 function F = f_local(einds,loc,grid,params)
250 %function F = f_local(einds,loc,grid,params)
252 % function to be evaluated in local coordinates, i.e.
253 % glob = [X, Y] with X and Y column vectors of global coordinates,
254 % producing a corresponding
255 % result-matrix in F (each row one result)
256 % such functions can also be used with 3-dim coordinates without
259 F = zeros(size(einds,2),dimrange);
261 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
264 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
265 function F = f_global(glob,params)
266 %function F = f_global(glob,params)
267 % function to be evaluated in gobal coordinates, i.e.
268 % glob = [X, Y] with X and Y column vectors of global coordinates,
269 % producing a corresponding
270 % result-matrix in F (each row one result)
271 % such functions can also be used with 3-dim coordinates without
274 F = zeros(size(glob,1),dimrange);
275 F(:,1) = sin(pi*glob(:,1));
276 F(:,2) = cos(pi*glob(:,2));
277 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A triangular conforming grid in two dimensions.
structure for the information of a ldg-function.
an ldg shape functions implementation
function demo_ldgfunc()
demo for showing ldgfunc capabilities