Add inner loop for multiple kacanov iterations

experiments/iterativeKacanov.m

@@ -7,6 +7,7 @@
 nEmax = 1e4;
 theta = 0.5;
 coeff = 1;
+lambda = 0.1;
 
 %% initialization
 mesh = Mesh.loadFromGeometry('unitsquare');
@@ -30,29 +31,42 @@ P = LoFeProlongation(fes);
 
 %% adaptive loop
 ell = 0;
-[nElements, estimator, iterationDiff, exactDiff] = deal(zeros(1000,1));
+[nElements, estimator, exactDiff, nIter] = deal(zeros(1000,1));
 meshSufficientlyFine = false;
 while ~meshSufficientlyFine
     ell = ell + 1;
+    freeDofs = getFreeDofs(fes);
     A = assemble(blf);
-    F = assemble(lf);
+    A = A(freeDofs,freeDofs);
     
     K = assemble(perturbation);
-    Aref = A + K;
+    Aref = A + K(freeDofs,freeDofs);
     Fref = assemble(lfRef);
+    uref.setFreeData(Aref \ Fref(freeDofs));
     
-    freeDofs = getFreeDofs(fes);
-    u.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
-    uref.setFreeData(Aref(freeDofs,freeDofs) \ Fref(freeDofs));
+    % inner iteration: sufficiently many kacanov iterations
+    n = 0;
+    solutionAccurateEnough = false;
+    while ~solutionAccurateEnough
+        n = n + 1;
+        F = assemble(lf);
+        u.setFreeData(A \ F(freeDofs));
+        eta2 = estimate(blf, lf, u);
+        
+        v = u.data(freeDofs) - uold.data(freeDofs);
+        iterationDiff = sqrt(v * A * v');
+        
+        uold.setData(u.data);
+        solutionAccurateEnough = (iterationDiff < lambda*sqrt(sum(eta2)));
+    end
     
-    eta2 = estimate(blf, lf, u);
     nElements(ell) = mesh.nElements;
     estimator(ell) = sqrt(sum(eta2));
-    v = u.data - uold.data;
-    iterationDiff(ell) = sqrt(v * A * v');
-    v = u.data - uref.data;
+    v = u.data(freeDofs) - uref.data(freeDofs);
     exactDiff(ell) = sqrt(v * A * v');
-    printLogMessage('number of elements: %d, estimator: %.2e', mesh.nElements, estimator(ell));
+    nIter(ell) = n;
+    printLogMessage('number of elements: %d, estimator: %.2e after %d iterations', ...
+        mesh.nElements, estimator(ell), n);
     
     meshSufficientlyFine = (mesh.nElements > nEmax);
     
@@ -70,15 +84,19 @@ end
 
 figure()
 idx = 1:ell;
+yyaxis left
+ylabel('error / estimator')
 loglog(nElements(idx), estimator(idx), '-o', 'LineWidth', 2)
 hold on
-loglog(nElements(idx), iterationDiff(idx), '-o', 'LineWidth', 2)
 loglog(nElements(idx), exactDiff(idx), '-o', 'LineWidth', 2)
 shift = estimator(1)*nElements(1);
 loglog(nElements(idx), shift*nElements(idx).^(-1/2), '--', 'LineWidth', 2, 'Color', 'k')
-legend('\eta_H', '|u^k_H - u^{k-1}_H|_{H^1}', '|u_H - u^{k}_H|_{H^1}', '#T_H^{-1/2}')
+yyaxis right
+semilogx(nElements(idx), nIter(idx), '-x', 'LineWidth', 2)
+ylim([0, max(nIter)+1])
+legend('\eta_H', '|u_H - u^{k}_H|_{H^1}', '#T_H^{-1/2}', 'iteration count')
+ylabel('number of iterations')
 xlabel('#T_H')
-ylabel('error / estimator')
 title('Kacanov iteration for non-symmetric / indefinite problem')