diff --git a/examples/convergenceRates.m b/paper-specific/convergenceRates.m
similarity index 89%
rename from examples/convergenceRates.m
rename to paper-specific/convergenceRates.m
index 2e200e793e389b7501548882e6a5b6b7bbb99103..796c76ac54bb81ec835b5bc939c10de1ec2d696a 100644
--- a/examples/convergenceRates.m
+++ b/paper-specific/convergenceRates.m
@@ -1,12 +1,13 @@
% ******************************************************************************
% Adaptive FEM algorithm with known solution and convergence rates for higher
% order finite elements.
+% This was used to generate data for Figure 6.
% ******************************************************************************
%% paramters
-nDofsMax = 1e4;
+nDofsMax = 1e7;
theta = 0.5;
-pmax = 8;
+pmax = 4;
%% initialization for every polynomial degree
[nElem, nDofs, errEst, h1Err] = deal(zeros(pmax, 1000));
@@ -60,9 +61,18 @@ for p = 1:pmax
end
end
-%% plot convergence rates
-plotData(nDofs, 'number of dofs', errEst, 'error estimator')
-plotData(nDofs, 'number of dofs', h1Err, 'H^1 error')
+%% plot convergence ratesofs, 'number of dofs', h1Err, 'H^1 error')
+% plotData(nDofs, 'number of dofs', errEst, 'error estimator')
+% plotData(nDofs, 'number of dofs', h1Err, 'H^1 error')
+
+for p = 1:pmax
+ idx = find(nElem(p,:) > 0);
+ fileID = fopen(['afemP', num2str(p), '.dat'], 'w');
+ fprintf(fileID, 'nElements,nDofs,estimator,H1Error\n');
+ fprintf(fileID, '%d,%d,%.3e,%.3e\n', ...
+ [nElem(p,idx); nDofs(p,idx); errEst(p,idx); h1Err(p,idx)]);
+ fclose(fileID);
+end
%% helper function for plotting
function plotData(xdata, xlab, ydata, ylab)
diff --git a/examples/goalOrientedAFEM.m b/paper-specific/goalOrientedAFEM.m
similarity index 81%
rename from examples/goalOrientedAFEM.m
rename to paper-specific/goalOrientedAFEM.m
index c270eba74963aad09d60c3e803885ae2a59e6af8..205e2c3223af4e2ed772e969bda5f539552117bc 100644
--- a/examples/goalOrientedAFEM.m
+++ b/paper-specific/goalOrientedAFEM.m
@@ -1,9 +1,10 @@
% ******************************************************************************
% Example of a GOAFEM algorithm taken from [Mommer, Stevenson; 2009].
+% This was used to generate data for Figure 8 (left).
% ******************************************************************************
%% paramters
-nDofsMax = 1e4;
+nDofsMax = 1e7;
theta = 0.5;
pmax = 3;
@@ -76,21 +77,29 @@ for p = 1:pmax
end
%% plot convergence rates
-pmax = size(nDofs, 1);
-figure()
+% figure()
+% for p = 1:pmax
+% idx = find(nDofs(p,:) > 0);
+% loglog(nDofs(p,idx), goalErrEst(p,idx), '-o', 'LineWidth', 2, 'DisplayName', ['p=',num2str(p)]);
+% hold on
+% end
+% for p = unique([1,pmax])
+% x = nDofs(p,idx) / nDofs(p,1);
+% loglog(nDofs(p,1)*x, goalErrEst(p,1)*x.^(-p), '--', 'LineWidth', 2, 'Color', 'k', 'DisplayName', ['\alpha = ', num2str(p)])
+% end
+% legend
+% xlabel('number of dofs')
+% ylabel('goal error estimator')
+% title('goal error estimator over number of dofs')
+
for p = 1:pmax
- idx = find(nDofs(p,:) > 0);
- loglog(nDofs(p,idx), goalErrEst(p,idx), '-o', 'LineWidth', 2, 'DisplayName', ['p=',num2str(p)]);
- hold on
-end
-for p = unique([1,pmax])
- x = nDofs(p,idx) / nDofs(p,1);
- loglog(nDofs(p,1)*x, goalErrEst(p,1)*x.^(-p), '--', 'LineWidth', 2, 'Color', 'k', 'DisplayName', ['\alpha = ', num2str(p)])
+ idx = find(nElem(p,:) > 0);
+ fileID = fopen(['goafemP', num2str(p), '.dat'], 'w');
+ fprintf(fileID, 'nElements,nDofs,goalErrorEstimate\n');
+ fprintf(fileID, '%d,%d,%.3e\n', ...
+ [nElem(p,idx); nDofs(p,idx); goalErrEst(p,idx)]);
+ fclose(fileID);
end
-legend
-xlabel('number of dofs')
-ylabel('goal error estimator')
-title(['goal error estimator over number of dofs'])
%% local function for residual a posteriori error estimation
% \eta(T)^2 = h_T^2 * || \Delta u ||_{L^2(T)}^2
diff --git a/examples/iterativeLinearization.m b/paper-specific/iterativeLinearization.m
similarity index 89%
rename from examples/iterativeLinearization.m
rename to paper-specific/iterativeLinearization.m
index 7e199187ec144f0b0ec32084bb59f5c9ea6072bb..c4f988eb2ff2e14cb3e12f8e8546c82f1d44c516 100644
--- a/examples/iterativeLinearization.m
+++ b/paper-specific/iterativeLinearization.m
@@ -1,16 +1,17 @@
% ******************************************************************************
% Example from Section 5.3 of [Heid, Praetorius, Wihler; 2021].
+% This was used to generate data for Figure 8 (right).
% ******************************************************************************
%% paramters
-nElemMax = 1e4;
+nElemMax = 1e7;
theta = 0.5;
lambda = 0.5;
deltaZ = 0.1;
linearizations = ["zarantonello", "kacanov", "newton"];
%% initialization for every linearization method
-[nElem, eta, E] = deal(zeros(length(linearizations), 1000));
+[nElem, eta, E, time] = deal(zeros(length(linearizations), 1000));
for k = 1:length(linearizations)
%% setup geometry & spaces
printLogMessage('*** starting %s iteration', linearizations(k))
@@ -27,6 +28,7 @@ for k = 1:length(linearizations)
u.setData(0);
%% adaptive algorithm
+ tic
i = 1;
Eold = 0;
qr = QuadratureRule.ofOrder(1);
@@ -51,6 +53,7 @@ for k = 1:length(linearizations)
if i > 1, Eold = E(k,i-1); end
E(k,i) = sum(integrateElement(eDensity, qr)) + u.data * F;
+ time(k,i) = toc;
%% stopping criterion of the iterative solver
if sqrt(E(k,i) - Eold) <= lambda * eta(k,i)
@@ -79,7 +82,17 @@ for k = 1:length(linearizations)
end
%% plot convergence rates
-plotData(nElem, 'number of elements', eta, 'error estimator', linearizations)
+% plotData(nElem, 'number of elements', eta, 'error estimator', linearizations)
+% plotData(time, 'computation time', eta, 'error estimator', linearizations)
+
+for k = 1:length(linearizations)
+ idx = find(nElem(k,:) > 0);
+ fileID = fopen(strcat('ailfem-', linearizations(k), '.dat'), 'w');
+ fprintf(fileID, 'nElements,estimator,energy,time\n');
+ fprintf(fileID, '%d,%.3e,%.3e,%.3e\n', ...
+ [nElem(k,idx); eta(k,idx); E(k,idx); time(k,idx)]);
+ fclose(fileID);
+end
%% helper function for plotting
function plotData(xdata, xlab, ydata, ylab, names)
diff --git a/paper-specific/listing1.m b/paper-specific/listing1.m
new file mode 100644
index 0000000000000000000000000000000000000000..9aaab0543bae0d40dbdf51e7cdb968d6246d85ff
--- /dev/null
+++ b/paper-specific/listing1.m
@@ -0,0 +1,30 @@
+% ******************************************************************************
+% Simple AFEM algorithm for the Poisson equation (Listing 1 in the
+% accompanying paper).
+% ******************************************************************************
+
+mesh = Mesh.loadFromGeometry('unitsquare');
+fes = FeSpace(mesh, LowestOrderH1Fe);
+u = FeFunction(fes);
+blf = BilinearForm(fes);
+blf.a = Constant(mesh, 1);
+lf = LinearForm(fes);
+lf.f = Constant(mesh, 1);
+
+while mesh.nElements < 1e6
+ A = assemble(blf);
+ F = assemble(lf);
+ freeDofs = getFreeDofs(fes);
+ u.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
+
+ hT = sqrt(getAffineTransformation(mesh).area);
+ qrEdge = QuadratureRule.ofOrder(1, '1D');
+ qrTri = QuadratureRule.ofOrder(1, '2D');
+ volumeRes = integrateElement(CompositeFunction(@(x) x.^2, lf.f), qrTri);
+ edgeRes = integrateNormalJump(Gradient(u), qrEdge, @(j) j.^2, {}, ':');
+ edgeRes(mesh.boundaries{:}) = 0;
+ eta2 = hT.^2.*volumeRes + hT.*sum(edgeRes(mesh.element2edges), Dim.Vector);
+
+ marked = markDoerflerSorting(eta2, 0.5);
+ mesh.refineLocally(marked, 'NVB');
+end
\ No newline at end of file
diff --git a/paper-specific/resolveComputationTime.m b/paper-specific/resolveComputationTime.m
new file mode 100644
index 0000000000000000000000000000000000000000..de971f158e2cf34883231f896ec5806d9c41c9eb
--- /dev/null
+++ b/paper-specific/resolveComputationTime.m
@@ -0,0 +1,133 @@
+% ******************************************************************************
+% Adaptive FEM algorithm with known solution and convergence rates.
+% This was used to generate data for Figure 7.
+% ******************************************************************************
+
+%% paramters
+nDofsMax = 1e6;
+theta = 0.5;
+p = 1;
+
+%% initialization
+[nElem, nDofs, errEst] = deal(zeros(1, 1000));
+time = zeros(6,1000);
+
+%% setup geometry & spaces
+printLogMessage('*** p = %d ***', p)
+mesh = Mesh.loadFromGeometry('Lshape');
+fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', 1);
+u = FeFunction(fes);
+uex = FeFunction(fes);
+
+%% set problem data for -\Delta u = 1 on L-shape
+blf = BilinearForm(fes);
+lf = LinearForm(fes);
+
+blf.a = Constant(mesh, 1);
+blf.qra = QuadratureRule.ofOrder(max(2*p-2, 1));
+
+lf.neumann = MeshFunction(mesh, @exactSolutionNeumannData);
+lf.qrNeumann = QuadratureRule.ofOrder(2*p, '1D');
+lf.bndNeumann = 2;
+
+%% adaptive loop
+i = 1;
+tic
+while 1
+ %% assemble & solve FEM system
+ A = assemble(blf);
+ time(1,i) = toc;
+ F = assemble(lf);
+ time(2,i) = toc;
+ free = getFreeDofs(fes);
+ u.setFreeData(A(free,free) \ F(free));
+ time(3,i) = toc;
+
+ %% estimate error and store data
+ eta2 = estimate(blf, lf, u);
+ errEst(i) = sqrt(sum(eta2));
+ time(4,i) = toc;
+ nDofs(i) = getDofs(fes).nDofs;
+ nElem(i) = mesh.nElements;
+ printLogMessage('number of dofs: %d, estimator: %.2e', nDofs(i), errEst(i));
+
+ %% refine mesh
+ marked = markDoerflerSorting(eta2, theta);
+ time(5,i) = toc;
+ mesh.refineLocally(marked, 'NVB');
+ time(6,i) = toc;
+
+ %% stoping criterion
+ if nDofs(i) > nDofsMax
+ break
+ end
+ i = i+1;
+end
+
+%% store fine grained computation times
+idx = find(nDofs > 0);
+time = time(:,idx);
+time = reshape([time(1,1); diff(time(:))], [], length(idx));
+fileID = fopen(['timingP', num2str(p), '.dat'], 'w');
+fprintf(fileID, 'nElements,nDofs,tAssembleA,tAssembleF,tSolve,tEstimate,tMark,tRefine,tCombined\n');
+fprintf(fileID, '%d,%d,%.3e,%.3e,%.3e,%.3e,%.3e,%.3e,%.3e\n', ...
+ [nElem(idx); nDofs(idx); time; sum(time, 1)]);
+fclose(fileID);
+
+%% local function for residual a posteriori error estimation
+function indicators = estimate(blf, lf, u)
+ p = blf.fes.finiteElement.order;
+ mesh = blf.fes.mesh;
+
+ % compute volume residual element-wise
+ if p == 1
+ volumeRes = 0;
+ else
+ f = CompositeFunction(@(D2u) (D2u(1,:,:) + D2u(4,:,:)).^2, Hessian(u));
+ qr = QuadratureRule.ofOrder(max(2*(p-2), 1));
+ volumeRes = integrateElement(f, qr);
+ end
+
+ % compute edge residual edge-wise
+ qr = QuadratureRule.ofOrder(p, '1D');
+ edgeRes = integrateNormalJump(Gradient(u), qr, ...
+ @(j) zeros(size(j)), {}, mesh.boundaries{1}, ... % no jump on dirichlet edges
+ @(j,phi) j-phi, {lf.neumann}, mesh.boundaries{2},...% neumann jump on neumann edges
+ @(j) j.^2, {}, ':'); % normal jump on inner edges
+
+ % combine the resdiuals suitably
+ hT = sqrt(getAffineTransformation(mesh).area);
+ indicators = hT.^2 .* volumeRes + ...
+ hT .* sum(edgeRes(mesh.element2edges), Dim.Vector);
+end
+
+%% exact solution for Lshape
+function y = exactSolution(x)
+ [phi, r] = cart2pol(x(1,:,:), x(2,:,:));
+ phi = phi + 2*pi*(phi < 0);
+ y = r.^(2/3) .* sin(2/3 * phi);
+end
+
+function y = exactSolutionNeumannData(x)
+ x1 = x(1,:,:);
+ x2 = x(2,:,:);
+
+ % determine boundary partsexactSolutionNeumannData
+ right = (x1 > 0) & (abs(x1) > abs(x2));
+ left = (x1 < 0) & (abs(x1) > abs(x2));
+ top = (x2 > 0) & (abs(x1) < abs(x2));
+ bottom = (x2 < 0) & (abs(x1) < abs(x2));
+
+ % compute d/dn uex
+ [phi, r] = cart2pol(x(1,:,:), x(2,:,:));
+ Cr = 2/3 * r.^(-4/3);
+ Cphi = 2/3 * (phi + 2*pi*(phi < 0));
+ dudx = Cr .* (x1.*sin(Cphi) - x2.*cos(Cphi));
+ dudy = Cr .* (x2.*sin(Cphi) + x1.*cos(Cphi));
+
+ y = zeros(size(x1));
+ y(right) = dudx(right);
+ y(left) =-dudx(left);
+ y(top) = dudy(top);
+ y(bottom) =-dudy(bottom);
+end