Clean up examples

examples/storingLevelOrientedData.m

-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;
     

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();