Add ZZ-estimator with Scott-Zhang as averaging operator

examples/zzEstimatorWithScottZhang.m

+% ******************************************************************************
+% Prototypical AFEM example with iterative solver for higher order finite
+% elements steered by ZZ-estimator with Scott-Zhang as averaging operator.
+% ******************************************************************************
+
+%% parameters
+nDofsMax = 1e5;
+theta = 0.5;
+pmax = 4;
+lamalg = 0.1;
+
+%% initialization for every polynomial degree
+[nElem, nDofs, errEst, iterations, time, directSolveTime] = deal(zeros(pmax, 1000));
+for p = 1:pmax
+    %% setup geometry & spaces
+    printLogMessage('*** p = %d (of %d) ***', p, pmax)
+    mesh = Mesh.loadFromGeometry('Lshape');
+    fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', ':');
+    u = FeFunction(fes);
+    u.setData(0);
+
+    %% set problem data for -\Delta u = 1 on L-shape
+    blf = BilinearForm();
+    lf = LinearForm();
+
+    blf.a = Constant(mesh, 1);
+    lf.f = Constant(mesh, 1);
+
+    %% set up solver and operator for nested iteration
+    % choose the iterative solver of choice: multigrid (with the variants
+    % lowOrderVcycle and highOrderVcycle), pcg (with jacobi and additive
+    % Schwarz preconditioner), and the cg solver
+    %[solver, P] = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+    %[solver, P] = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+    [solver, P] = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
+    solver.tol = 1e-6;
+    solver.maxIter = 100;
+
+    %% adaptive loop
+    ell = 0;
+    meshSufficientlyFine = false;
+    while ~meshSufficientlyFine
+        %% setup
+        tic;
+        ell = ell + 1;
+        A = assemble(blf, fes);
+        rhs = assemble(lf, fes);
+
+        freeDofs = getFreeDofs(fes);
+        A = A(freeDofs,freeDofs);
+        solver.setupSystemMatrix(A);
+        solver.setupRhs(rhs(freeDofs), u.data(freeDofs)');
+
+        % exact solution as reference
+        tic;
+        x_ex = A \ rhs(freeDofs);
+        directSolveTime(p, ell) = toc;
+        
+        %% iterative solve
+        % solver.solve() implements the algebraic solver loop and
+        % terminates automatically once the wanted accuracy is reached
+        % Another option is to manually implement a loop and implement one
+        % solver step with solver.step() and stop by setting a stopping
+        % criterion
+        solverIsConverged = false;
+        while ~solverIsConverged
+            solver.step();
+            u.setFreeData(solver.x)
+            mu2 = estimateZZ(blf, lf, u);
+            estimator = sqrt(sum(mu2));
+            solver.tol = lamalg*estimator;
+            solverIsConverged = solver.isConverged;
+        end
+        
+
+        %% estimate error and store data
+        nDofs(p,ell) = getDofs(fes).nDofs;
+        errEst(p,ell) = estimator;
+        iterations(p, ell) = solver.iterationCount();
+        nElem(p,ell) = mesh.nElements;
+        printLogMessage('number of dofs: %d, estimator: %.2e after %d iterations', ...
+            nDofs(p,ell), errEst(p,ell), solver.iterationCount);
+
+        %% stoping criterion
+        meshSufficientlyFine = (nDofs(p,ell) > nDofsMax);
+
+        %% refine mesh
+        if ~meshSufficientlyFine
+            marked = markDoerflerSorting(mu2, theta);
+            mesh.refineLocally(marked, 'NVB');
+            u.setData(prolongate(P, u));
+        end
+        time(p,ell) = toc;
+    end
+end
+
+figure()
+for p = 1:pmax
+    idx = find(nDofs(p,:) > 0);
+    complexity = [cumsum(nDofs(p, idx) .* iterations(p, idx))];
+    loglog(complexity, errEst(p,idx), '-o', 'LineWidth', 2, 'DisplayName', ['p=',num2str(p)]);
+    hold on
+end
+for p = unique([1,pmax])
+    x = complexity / complexity(1);
+    loglog(complexity(1)*x, errEst(p,1)*x.^(-p/2), '--', 'LineWidth', 2, 'Color', 'k', 'DisplayName', ['\alpha = ', num2str(p), '/2'])
+end
+legend
+xlabel('computational cost')
+ylabel('\eta')
+title('error estimator over number of total computational cost')
+
+
+%% local function for recovery based a posteriori error estimation
+
+function indicators = estimateZZ(~, ~, u)
+
+p = u.fes.finiteElement.order;
+qr = QuadratureRule.ofOrder(2*p);
+
+% Remark: Scott-Zhang operator for piecewise polynomials of degree p corresponds to nodal interpolation
+
+% nodal interpolation of first component of gradient of u
+gGradU1 = FeFunction(u.fes);
+gGradU1.setData(nodalInterpolation(CompositeFunction(@(grad) grad(1,:,:), Gradient(u)), u.fes));
+
+% nodal interpolation of second component of gradient of u
+gGradU2 = FeFunction(u.fes);
+gGradU2.setData(nodalInterpolation(CompositeFunction(@(grad) grad(2,:,:), Gradient(u)), u.fes));
+
+% compute recovery term
+recoveryTermHandle = @(g1, g2, grad) (g1 - grad(1,:,:)).^2 + (g2 - grad(2,:,:)).^2;
+recoveryTerm = CompositeFunction(recoveryTermHandle, gGradU1, gGradU2, Gradient(u));
+indicators = integrateElement(recoveryTerm, qr);
+
+end