diff --git a/examples/storingLevelOrientedData.m b/examples/storingLevelOrientedData.m
index 447203f12dc161c7a543f44bbd36dde0d59d8a00..5aa2c314dfaa28757354e995ce99c58f1709efc6 100644
--- a/examples/storingLevelOrientedData.m
+++ b/examples/storingLevelOrientedData.m
@@ -1,10 +1,10 @@
-function leveldata = storingLevelOrientedData(doPlots)
-%%MAINLEVELDATA Example of an adaptive FEM algorithm with storing and
-%displaying level-oriented dataa
+% Example of an adaptive FEM algorithm with storing and displaying
+% level-oriented data
- %% proceed optional input
- if nargin < 1
- doPlots = true;
+function leveldata = storingLevelOrientedData(doPlots, doStore)
+ arguments
+ doPlots (1,1) logical = true
+ doStore (1,1) logical = false
end
%% paramters
@@ -15,7 +15,7 @@ function leveldata = storingLevelOrientedData(doPlots)
mesh = Mesh.loadFromGeometry('Lshape');
fes = FeSpace(mesh, LowestOrderH1Fe());
u = FeFunction(fes);
-
+
%% initialize output data structure
pathToStorage = 'results';
leveldata = LevelData(pathToStorage);
@@ -23,69 +23,68 @@ function leveldata = storingLevelOrientedData(doPlots)
leveldata.metaData('domain') = 'Lshape';
leveldata.metaData('method') = 'S1';
leveldata.metaData('identifier') = 'example';
-
+
%% problem data
blf = BilinearForm();
lf = LinearForm();
- blf.a = Constant(mesh, [1;0;0;1]);
- blf.b = Constant(mesh, [0;0]);
- blf.c = Constant(mesh, 0);
+ blf.a = Constant(mesh, 1);
lf.f = Constant(mesh, 1);
- lf.fvec = Constant(mesh, [0;0]);
% reference value for energy norm of exact solution
normExactSquared = 0.2140750232;
-
+
%% adaptive loop
meshSufficientlyFine = false;
while ~meshSufficientlyFine
%% compute data for rhs coefficient and assemble forms
A = assemble(blf, fes);
F = assemble(lf, fes);
-
+
%% solve FEM system
freeDofs = getFreeDofs(fes);
tic;
u.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
runtimeSolve = toc;
-
+
%% approximate error
err = sqrt(normExactSquared - u.data * A * u.data');
-
+
%% compute refinement indicators
tic;
eta2 = estimate(blf, lf, u);
eta = sqrt(sum(eta2));
runtimeEstimate = toc;
-
+
%% compute efficiency index
eff = err / eta;
-
+
%% storing error quantities and general data
leveldata.append('ndof', uint32(length(freeDofs)), ...
- 'eta', eta, ...
- 'err', err, ...
- 'coordinates', mesh.coordinates, ...
- 'elements', mesh.elements);
-
+ 'eta', eta, ...
+ 'err', err, ...
+ 'coordinates', mesh.coordinates, ...
+ 'elements', mesh.elements);
+
%% storing results of timing
leveldata.setTime(leveldata.nLevel, ...
- 'runtimeSolve', runtimeSolve, ...
- 'runtimeEstimate', runtimeEstimate);
-
+ 'runtimeSolve', runtimeSolve, ...
+ 'runtimeEstimate', runtimeEstimate);
+
%% storing absolute values
leveldata.setAbsolute(leveldata.nLevel, 'eff', eff);
-
+
%% print information on current level to command line
leveldata.printLevel();
-
- %% save intermediate results to file
+
+ %% save intermediate results to file
% (prevents data loss in case of crashes)
- leveldata.saveToFile();
-
+ if doStore
+ leveldata.saveToFile();
+ end
+
%% stoping criterion
meshSufficientlyFine = (mesh.nElements > nEmax);
-
+
%% refine mesh
if ~meshSufficientlyFine
marked = markDoerflerBinning(eta2, theta);
@@ -95,8 +94,8 @@ function leveldata = storingLevelOrientedData(doPlots)
%% post-processing variables
leveldata.setTime(':', 'runtimeTotal', ...
- sum(leveldata.get(':', 'runtimeSolve', 'runtimeEstimate'), 2));
-
+ sum(leveldata.get(':', 'runtimeSolve', 'runtimeEstimate'), 2));
+
%% plot results
if doPlots
% plot error quantities (in general, converging to zeros)
@@ -110,54 +109,41 @@ function leveldata = storingLevelOrientedData(doPlots)
% plot absolute values in semi-logarithmic plot
figure();
leveldata.plotAbsolute('ndof');
-
- % plot of all triangulations
- % figure();
- % leveldata.plotTriangulation();
-
- %% export error plot to file
+ end
+
+ %% additional visualization/storage capabilities
+ if doStore
+ % export error plot to file
leveldata.plotToFile('ndof');
- %% export refinement video
- % (generates video of about 110 MB)
- %leveldata.plotTriangulationToFile();
- end
+ % save all values to comma-separated table
+ leveldata.saveToTable();
- %% save all values to comma-separated table
- leveldata.saveToTable();
+ % plot of all triangulations one after another
+ figure();
+ leveldata.plotTriangulation();
+
+ % export refinement video (generates video of about 110 MB)
+ leveldata.plotTriangulationToFile();
+ end
end
-
%% local function for residual a posteriori error estimation
-% This function shows how to build a residual error estimator with the provided
-% tools. Not hiding this in its own class has the advantage that estimation is
-% very flexible with respect to the involved terms.
-% \eta(T)^2 = h_T^2 * || -div(a*Du) + b*Du + c*u - f + div(fvec) ||_{L^2(T)}^2
-% + h_T * || [ a*Du - fvec ] ||_{L^2(E)}^2
-function indicators = estimate(blf, lf, u)
+% \eta(T)^2 = h_T^2 * || \Delta u + f ||_{L^2(T)}^2 + h_T * || [Du] ||_{L^2(E)}^2
+function indicators = estimate(~, ~, u)
fes = u.fes;
mesh = fes.mesh;
%% compute volume residual element-wise
- % Here, one can optimize the estimator for the given situation. In the case
- % p=1 with constant diffusion, the diffusion term vanishes in the residual.
- % Also div(fvec) vanishes if fvec is element-wise constant.
- f = CompositeFunction(...
- @(b, c, f, u, Du) (vectorProduct(b, Du) + c .* u - f).^2, ...
- blf.b, blf.c, lf.f, u, Gradient(u));
- qrTri = QuadratureRule.ofOrder(4);
+ f = CompositeFunction(@(Du) vectorProduct(Du, Du), Gradient(u));
+ qrTri = QuadratureRule.ofOrder(1);
volumeRes = integrateElement(f, qrTri);
%% compute edge residual edge-wise
- % Here, none of the terms vanish, so they must be computed. This can be done
- % with a lower order than the element-residuals.
- f = CompositeFunction(...
- @(a, fvec, Du) vectorProduct(a, Du, [2,2], [2,1]) - fvec, ...
- blf.a, lf.fvec, Gradient(u));
qrEdge = QuadratureRule.ofOrder(1, '1D');
- edgeRes = integrateNormalJump(f, qrEdge, @(j) j.^2, {}, ':');
+ edgeRes = integrateNormalJump(Gradient(u), qrEdge, @(j) j.^2, {}, ':');
dirichlet = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
edgeRes(dirichlet) = 0;
diff --git a/examples/timingWithLevelDataCollection.m b/examples/timingWithLevelDataCollection.m
index 95dcf89c9adb57d79f786633d945aab1a3dfb5de..ad36a75950a1d1099ad1442446b595734281ada8 100644
--- a/examples/timingWithLevelDataCollection.m
+++ b/examples/timingWithLevelDataCollection.m
@@ -8,7 +8,7 @@ nRuns = 10;
%% run time measurements
leveldatacollection = ...
- TimeIt('debugTiming', nRuns, 'storingLevelOrientedData', false);
+ TimeIt('debugTiming', nRuns, 'storingLevelOrientedData', false, false);
%% print statistical analysis of timing results
leveldatacollection.printStatistics();