diff --git a/.gitignore b/.gitignore
index 6e36dd534a1f46e0ea94f2e85cd478a7505a2f21..8ab9f475ea2b6aa42bfa029442d30929a35dccd2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,14 @@ slprj/
# Simulink autosave extension
.autosave
+**/.DS_Store
+
+# Ctags file
+tags
+
+# Vim files
+*.swp
+*.swo
+
+# Output files
+*.csv
diff --git a/README.md b/README.md
index ae98e01ca692969c03f285778519350260a8d83a..0e8f478f2b23fff53d5c32be3b8ca0b6fd485627 100644
--- a/README.md
+++ b/README.md
@@ -29,7 +29,7 @@ Advantages:
To get started, run `setup.m` in the root folder. This adds everything to the
path and compiles .mex files if necessary.
-Note that, at least, Matlab version R2020b is required.
+Note that, at least, Matlab version R2022b is required.
## First steps
diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..3ef04b1eb3a6f50650b74ff9ad1cf3380613f7b9
--- /dev/null
+++ b/examples/algebraicSolver.m
@@ -0,0 +1,143 @@
+% ******************************************************************************
+% Prototypical AFEM example with iterative solver for higher order finite
+% elements.
+% ******************************************************************************
+
+%% parameters
+nDofsMax = 1e4;
+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);
+ blf.qra = QuadratureRule.ofOrder(max(2*p-2, 1));
+
+ lf.f = Constant(mesh, 1);
+ lf.qrf = QuadratureRule.ofOrder(2*p);
+
+ %% 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, '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)
+ eta2 = estimate(blf, lf, u);
+ estimator = sqrt(sum(eta2));
+ 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(eta2, 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 residual a posteriori error estimation
+% \eta(T)^2 = h_T^2 * || \Delta u ||_{L^2(T)}^2
+% + h_T * || [[(Du - fvec) * n]] ||_{L^2(E) \cap \Omega}^2
+function indicators = estimate(blf, lf, u)
+p = u.fes.finiteElement.order;
+mesh = u.fes.mesh;
+trafo = getAffineTransformation(mesh);
+
+% compute volume residual element-wise
+% For p=1, the diffusion term vanishes in the residual.
+if p == 1
+ f = CompositeFunction(@(f) f.^2, lf.f);
+else
+ f = CompositeFunction(@(v,f) (v(1,:,:) + v(4,:,:) + f).^2, Hessian(u), lf.f);
+end
+qr = QuadratureRule.ofOrder(2*p);
+volumeResidual = integrateElement(f, qr);
+
+% compute edge residual edge-wise
+qr = QuadratureRule.ofOrder(max(p-1, 1), '1D');
+edgeResidual = integrateNormalJump(Gradient(u), qr, ...
+ @(j) zeros(size(j)), {}, vertcat(mesh.boundaries{:}), ... % no jump on dirichlet edges
+ @(j) j.^2, {}, ':'); % normal jump on inner edges
+
+% combine the residuals suitably
+hT = sqrt(trafo.area);
+indicators = hT.^2 .* volumeResidual + ...
+ hT .* sum(edgeResidual(mesh.element2edges), Dim.Vector);
+end
diff --git a/examples/goafemIterativeSolver.m b/examples/goafemIterativeSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..50a7df3835d03592ed0f894e33081a62b8935177
--- /dev/null
+++ b/examples/goafemIterativeSolver.m
@@ -0,0 +1,150 @@
+% ******************************************************************************
+% Example of a GOAFEM algorithm with iterative solver.
+% ******************************************************************************
+
+%% paramters
+nDofsMax = 1e4;
+theta = 0.5;
+pmax = 4;
+lamalg = 0.1;
+
+%% setup geometry & spaces
+[nElem, nDofs, nIterPrimal, nIterDual, goalErrEst] = deal(zeros(1, 1000));
+
+for p = 1:pmax
+ printLogMessage('*** GOAFEM with p = %d and iterative solver ***', p)
+ mesh = Mesh.loadFromGeometry('unitsquare');
+ mesh.refineUniform(1, 'RGB');
+ fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', ':');
+ u = FeFunction(fes);
+ z = FeFunction(fes);
+ u.setData(0);
+ z.setData(0);
+
+ %% set problem data
+ % primal: -\Delta u = "lorentzian peak 1"
+ % dual: -\Delta z = "lorentzian peak 2"
+ blf = BilinearForm();
+ blf.a = Constant(mesh, 1);
+ blf.qra = QuadratureRule.ofOrder(max(2*p-2, 1));
+
+ lfF = LinearForm();
+ lfF.f = MeshFunction(mesh, @(x) lorentzian(x, [0.7;0.7], 1e-1));
+ lfF.qrf = QuadratureRule.ofOrder(2*p);
+
+ lfG = LinearForm();
+ lfG.f = MeshFunction(mesh, @(x) lorentzian(x, [0.2;0.3], 1e-2));
+ lfG.qrf = QuadratureRule.ofOrder(2*p);
+
+ %% set up solver and lifting operator for nested iteration
+ [solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
+ solver.tol = 1e-6;
+ solver.maxIter = 1000;
+
+ %% adaptive loop
+ ell = 0;
+ meshSufficientlyFine = false;
+ while ~meshSufficientlyFine
+ ell = ell + 1;
+
+ %% assemble & solve FEM system iteratively
+ freeDofs = getFreeDofs(fes);
+ A = assemble(blf, fes);
+ A = A(freeDofs, freeDofs);
+ rhs = [assemble(lfF, fes), assemble(lfG, fes)];
+ rhs = rhs(freeDofs,:);
+ uz0 = [u.data', z.data'];
+ solver.setupSystemMatrix(A);
+ solver.setupRhs(rhs, uz0(freeDofs,:));
+
+ solverIsConverged = false;
+ while ~solverIsConverged
+ solver.step();
+ u.setFreeData(solver.x(:,1));
+ z.setFreeData(solver.x(:,2));
+ eta2 = estimate(blf, lfF, u);
+ zeta2 = estimate(blf, lfG, z);
+ primalEstimator = sqrt(sum(eta2));
+ dualEstimator = sqrt(sum(zeta2));
+ solver.tol = lamalg*[primalEstimator, dualEstimator];
+ solverIsConverged = solver.isConverged;
+ end
+ u.setFreeData(solver.x(:,1));
+ z.setFreeData(solver.x(:,2));
+
+ %% estimate error and store data
+ eta2 = estimate(blf, lfF, u);
+ zeta2 = estimate(blf, lfG, z);
+ nDofs(ell) = getDofs(fes).nDofs;
+ nElem(ell) = mesh.nElements;
+ nIterPrimal(ell) = solver.iterationCount(1);
+ nIterDual(ell) = solver.iterationCount(2);
+ goalErrEst(ell) = primalEstimator*dualEstimator;
+ printLogMessage('number of dofs: %d, iterations (primal/dual): %d / %d, estimator: %.2e', ...
+ nDofs(ell), nIterPrimal(ell), nIterDual(ell), goalErrEst(ell));
+
+ %% stoping criterion
+ meshSufficientlyFine = (nDofs(ell) > nDofsMax);
+
+ %% refine mesh and transfer solutions to finer mesh for nested iteration
+ if ~meshSufficientlyFine
+ marked = markGoafemMS(eta2, zeta2, theta);
+ mesh.refineLocally(marked, 'NVB');
+ u.setData(prolongate(P, u));
+ z.setData(prolongate(P, z));
+ end
+ end
+end
+
+%% plot convergence rates
+figure()
+idx = 1:ell;
+x = cumsum(nDofs(idx).*max(nIterPrimal(idx), nIterDual(idx)));
+yyaxis left
+loglog(x, goalErrEst(idx), '-x', 'LineWidth', 2, 'MarkerSize', 8);
+hold on
+loglog(x, goalErrEst(1)*(x/x(1)).^(-p), '--', 'LineWidth', 2, 'Color', 'k')
+ylabel('goal error estimator of final iterates')
+yyaxis right
+stem(x, nIterPrimal(idx), 's', 'filled')
+stem(x, nIterDual(idx), 'o', 'filled')
+ylabel('number of iterations')
+legend('\eta_l \zeta_l', ['\alpha=', num2str(p)], 'primal iterations', 'dual iteration')
+xlabel('total work')
+title(['GOAFEM p=', num2str(p)])
+
+%% local function for residual a posteriori error estimation
+% \eta(T)^2 = h_T^2 * || \Delta u ||_{L^2(T)}^2
+% + h_T * || [[(Du - fvec) * n]] ||_{L^2(E) \cap \Omega}^2
+function indicators = estimate(blf, lf, u)
+ p = u.fes.finiteElement.order;
+ mesh = u.fes.mesh;
+ trafo = getAffineTransformation(mesh);
+
+ % compute volume residual element-wise
+ % For p=1, the diffusion term vanishes in the residual.
+ if p == 1
+ f = lf.f;
+ else
+ f = CompositeFunction(@(v,f) (v(1,:,:) + v(4,:,:) + f).^2, Hessian(u), lf.f);
+ end
+ qr = QuadratureRule.ofOrder(2*p);
+ volumeResidual = integrateElement(f, qr);
+
+ % compute edge residual edge-wise
+ qr = QuadratureRule.ofOrder(max(p-1, 1), '1D');
+ edgeResidual = integrateNormalJump(Gradient(u), qr, ...
+ @(j) zeros(size(j)), {}, vertcat(mesh.boundaries{:}), ... % no jump on dirichlet edges
+ @(j) j.^2, {}, ':'); % normal jump on inner edges
+
+ % combine the resdiuals suitably
+ hT = sqrt(trafo.area);
+ indicators = hT.^2 .* volumeResidual + ...
+ hT .* sum(edgeResidual(mesh.element2edges), Dim.Vector);
+end
+
+% Lorentz-peak 1/(eps+||x-x0||^2)
+function y = lorentzian(x, x0, eps)
+ dx = x - x0;
+ y = 1 ./ (eps + vectorProduct(dx,dx));
+end
diff --git a/examples/goalOrientedAFEM.m b/examples/goalOrientedAFEM.m
index 6207bb6c78a6f13afc91c3fddf1c577c6e73f45f..69297791ea0e34e073a1890ebafb6db2d40e3140 100644
--- a/examples/goalOrientedAFEM.m
+++ b/examples/goalOrientedAFEM.m
@@ -39,7 +39,7 @@ for p = 1:pmax
lfG.qrfvec = QuadratureRule.ofOrder(max(p-1, 1));
%% set up lifting operators for rhs FEM-data
- P = LoFeProlongation(ncFes);
+ P = LoMeshProlongation(ncFes);
%% adaptive loop
ell = 1;
diff --git a/examples/iterativeLinearization.m b/examples/iterativeLinearization.m
index d37f9228b165a009ee53c9584bc827148038bbad..e4bf7d3b7c2dfe85e1f7485ab22b24c5f5aa97c6 100644
--- a/examples/iterativeLinearization.m
+++ b/examples/iterativeLinearization.m
@@ -27,7 +27,7 @@ for k = 1:length(linearizations)
u.setData(0);
%% nested iteration
- P = LoFeProlongation(fes);
+ P = LoMeshProlongation(fes);
%% adaptive algorithm
ell = 0;
diff --git a/examples/storingLevelOrientedData.m b/examples/storingLevelOrientedData.m
new file mode 100644
index 0000000000000000000000000000000000000000..74888820894fc9ebda98813c71343a9d20efdd8d
--- /dev/null
+++ b/examples/storingLevelOrientedData.m
@@ -0,0 +1,154 @@
+% Example of an adaptive FEM algorithm with storing and displaying
+% level-oriented data
+
+function leveldata = storingLevelOrientedData(doPlots, doStore)
+ arguments
+ doPlots (1,1) logical = true
+ doStore (1,1) logical = false
+ end
+
+ %% paramters
+ nEmax = 1e4;
+ theta = 0.5;
+
+ %% initialization
+ mesh = Mesh.loadFromGeometry('Lshape');
+ fes = FeSpace(mesh, LowestOrderH1Fe());
+ u = FeFunction(fes);
+
+ %% initialize output data structure
+ pathToStorage = 'results';
+ leveldata = LevelData(pathToStorage);
+ leveldata.metaData('problem') = 'Poisson';
+ leveldata.metaData('domain') = 'Lshape';
+ leveldata.metaData('method') = 'S1';
+ leveldata.metaData('identifier') = 'example';
+
+ %% problem data
+ blf = BilinearForm();
+ lf = LinearForm();
+ blf.a = Constant(mesh, 1);
+ lf.f = Constant(mesh, 1);
+
+ % 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);
+
+ %% storing results of timing
+ leveldata.setTime(leveldata.nLevel, ...
+ '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
+ % (prevents data loss in case of crashes)
+ if doStore
+ leveldata.saveToFile();
+ end
+
+ %% stoping criterion
+ meshSufficientlyFine = (mesh.nElements > nEmax);
+
+ %% refine mesh
+ if ~meshSufficientlyFine
+ marked = markDoerflerBinning(eta2, theta);
+ mesh.refineLocally(marked, 'NVB');
+ end
+ end
+
+ %% post-processing variables
+ leveldata.setTime(':', 'runtimeTotal', ...
+ sum(leveldata.get(':', 'runtimeSolve', 'runtimeEstimate'), 2));
+
+ %% plot results
+ if doPlots
+ % plot error quantities (in general, converging to zeros)
+ figure();
+ leveldata.plot('ndof');
+
+ % plot runtimes (in general, diverging to infinity)
+ figure();
+ leveldata.plotTime('ndof');
+
+ % plot absolute values in semi-logarithmic plot
+ figure();
+ leveldata.plotAbsolute('ndof');
+ end
+
+ %% additional visualization/storage capabilities
+ if doStore
+ % export error plot to file
+ leveldata.plotToFile('ndof');
+
+ % 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 2MB)
+ leveldata.plotTriangulationToFile();
+ end
+end
+
+
+
+%% local function for residual a posteriori error estimation
+% \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
+ f = CompositeFunction(@(Du) vectorProduct(Du, Du), Gradient(u));
+ qrTri = QuadratureRule.ofOrder(1);
+ volumeRes = integrateElement(f, qrTri);
+
+ %% compute edge residual edge-wise
+ qrEdge = QuadratureRule.ofOrder(1, '1D');
+ edgeRes = integrateNormalJump(Gradient(u), qrEdge, @(j) j.^2, {}, ':');
+ dirichlet = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
+ edgeRes(dirichlet) = 0;
+
+ %% combine the resdiuals suitably
+ hT = sqrt(getAffineTransformation(mesh).area);
+ indicators = hT.^2 .* volumeRes + ...
+ hT .* sum(edgeRes(mesh.element2edges), Dim.Vector);
+end
diff --git a/examples/timingWithLevelDataCollection.m b/examples/timingWithLevelDataCollection.m
new file mode 100644
index 0000000000000000000000000000000000000000..720e65e52a421384ed721ef1889de8db12a39fe0
--- /dev/null
+++ b/examples/timingWithLevelDataCollection.m
@@ -0,0 +1,22 @@
+% ******************************************************************************
+% Example of timing with LevelDataCollection
+% ******************************************************************************
+
+
+%% parameters
+nRuns = 10;
+
+%% run time measurements
+storeResults = false;
+leveldatacollection = ...
+ TimeIt('debugTiming', nRuns, storeResults, 'storingLevelOrientedData', false, false);
+
+%% print statistical analysis of timing results
+leveldatacollection.printStatistics();
+
+%% plot results
+figure();
+leveldatacollection.plotStatistics('ndof');
+
+%% save results to comma-separated file
+% leveldatacollection.saveToTable();
diff --git a/experiments/p_laplace.m b/experiments/p_laplace.m
new file mode 100644
index 0000000000000000000000000000000000000000..9bee4bdbb032684aa188423747fe374f545fd6f1
--- /dev/null
+++ b/experiments/p_laplace.m
@@ -0,0 +1,172 @@
+%% paramters
+nElemMax = 1e5;
+theta = 0.5;
+lambda = 0.1;
+deltaZ = 0.1;
+linearization = 'kacanov';
+geometry = 'Lshape';
+p = 3;
+
+%% initialisation
+[nElem, eta, E] = deal(zeros(1,1000));
+
+%% setup geometry & spaces
+mesh = Mesh.loadFromGeometry(geometry);
+mesh.refineUniform(2);
+fes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', ':');
+u = FeFunction(fes);
+
+%% set problem data given in [Heid, Praetorius, Wihler; 2021]
+Du = Gradient(u);
+g = Constant(mesh, 1);
+[blf, lf] = setProblemData(mesh, fes, Du, g, linearization, p);
+eDensity = CompositeFunction(@(t) muIntegral(vectorProduct(t, t),p), Du);
+u.setFreeData(rand(1,length(getFreeDofs(u.fes))));
+
+%% nested iteration
+P = LoMeshProlongation(fes);
+
+%% adaptive algorithm
+i = 1;
+Eold = 0;
+qr = QuadratureRule.ofOrder(1); % tweak?
+while 1
+ % the volume residual for lowest order is independent of the approximate
+ % solution, so it can be computed once
+ vol = estimateVolume(g);
+
+ %% inner iteration to solve the nonlinear problem
+ nIterations = 1;
+ while true
+ %% solve linearized problem
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
+ updateDataU(u, deltaZ, A, F, linearization);
+
+ %% compute remainder of error estimator and energy
+ edge = estimateEdge(u, mesh, p);
+ eta2 = combineEstimators(vol, edge, mesh);
+ eta(i) = sqrt(sum(eta2));
+ nElem(i) = mesh.nElements;
+
+ if i > 1, Eold = E(i-1); end
+ E(i) = sum(integrateElement(eDensity, qr)) + u.data * F;
+
+ %% stopping criterion of the iterative solver
+ if sqrt(E(i) - Eold) <= lambda * eta(i) % norm(e) <= fehler
+ break
+ end
+ nIterations = nIterations + 1;
+ i = i+1;
+ end
+
+ %% stoping criterion for the AFEM algorithm
+ printLogMessage('number of elements: %d, iterations: %d, estimator: %.2e', ...
+ nElem(i), nIterations, eta(i));
+ if nElem(i) > nElemMax
+ break
+ end
+
+ %% refine mesh
+ marked = markDoerflerSorting(eta2, theta);
+ mesh.refineLocally(marked, 'NVB');
+ u.setData(prolongate(P, u));
+ i = i+1;
+end
+
+%% plot convergence rates
+plotData(nElem, 'number of elements', eta, 'error estimator', linearization,p)
+
+%% helper function for plotting
+function plotData(xdata, xlab, ydata, ylab, name,p)
+ figure()
+ for k = 1:size(xdata,1)
+ idx = find(xdata(k,:) > 0);
+ loglog(xdata(k,idx), ydata(k,idx), '-o', 'LineWidth', 2, 'DisplayName', strcat(name,', p=', num2str(p)));
+ hold on
+ end
+ x = xdata(1,idx) / xdata(1,1);
+ loglog(xdata(1,1)*x, ydata(1,1)*x.^(-1/2), '--', 'LineWidth', 2, 'Color', 'k', 'DisplayName', '\alpha = 1/2')
+ legend
+ xlabel(xlab)
+ ylabel(ylab)
+ title([ylab, ' over ', xlab])
+end
+
+%% local function for problem data
+function [blf, lf] = setProblemData(mesh, fes, Du, g, linearization, p)
+ blf = BilinearForm();
+ lf = LinearForm();
+
+ switch linearization
+ case "zarantonello"
+ blf.a = Constant(mesh, 1);
+ lf.f = g;
+ lf.fvec = CompositeFunction(@(t) -mu(vectorProduct(t, t),p) .* t, Du);
+ case "kacanov"
+ blf.a = CompositeFunction(@(t) mu(vectorProduct(t, t),p), Du);
+ lf.f = g;
+ lf.fvec = [];
+ case "newton"
+ blf.a = CompositeFunction(@(t) mu(vectorProduct(t, t),p) .* [1;0;0;1] ...
+ + 2*muDerivative(vectorProduct(t, t),p).*vectorProduct(t, t, [2,1], [2,1]'), Du);
+ lf.f = g;
+ lf.fvec = CompositeFunction(@(t) -mu(vectorProduct(t, t),p) .* t, Du);
+ end
+
+ [blf.qra, lf.qrf, lf.qrfvec] = deal(QuadratureRule.ofOrder(1));
+end
+
+function updateDataU(u, deltaZ, A, F, linearization)
+ freeDofs = getFreeDofs(u.fes);
+ v = FeFunction(u.fes);
+ switch linearization
+ case "zarantonello"
+ v.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
+ u.setData(u.data + deltaZ*v.data);
+ case "kacanov"
+ %u.setData(A \ F);
+ u.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
+ case "newton"
+ v.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
+ u.setData(u.data + v.data);
+ end
+end
+
+%% local function for residual a posteriori error estimation
+% \eta(T)^2 = h_T^2 * || g ||_{L^2(T)}^2
+% + h_T * || [[ mu(|Du|^2)*Du * n]] ||_{L^2(E) \cap \Omega}^2
+function vol = estimateVolume(g)
+ qr = QuadratureRule.ofOrder(1);
+ vol = integrateElement(g, qr);
+end
+
+function edge = estimateEdge(u, mesh, p)
+ qr = QuadratureRule.ofOrder(1, '1D');
+ dirichlet = vertcat(mesh.boundaries{:});
+
+ f = CompositeFunction(@(t) mu(vectorProduct(t, t),p).*t, Gradient(u));
+ edge = integrateNormalJump(f, qr, ...
+ @(j) zeros(size(j)), {}, dirichlet, ... % no jump on dirichlet edges
+ @(j) j.^2, {}, ':'); % normal jump on inner edges
+
+end
+
+function indicators = combineEstimators(vol, edge, mesh)
+ trafo = getAffineTransformation(mesh);
+ hT = sqrt(trafo.area);
+ indicators = hT.^2.*vol + hT.*sum(edge(mesh.element2edges), Dim.Vector);
+end
+
+%% nonlinearities
+function y = mu(t,p)
+ y = t.^(p/2 - 1);
+end
+
+function y = muIntegral(t,p)
+ y = p/4 * t.^(p/2);
+end
+
+function y = muDerivative(t,p)
+ y = (p/2 - 1) * t.^(p/2 - 2);
+end
diff --git a/lib/assembly/@BilinearForm/BilinearForm.m b/lib/assembly/@BilinearForm/BilinearForm.m
index fc76b96c640553918080f09d79386e6a1e0bc42e..1221dd7053b48b1cea9af6cce316a872aa0add83 100644
--- a/lib/assembly/@BilinearForm/BilinearForm.m
+++ b/lib/assembly/@BilinearForm/BilinearForm.m
@@ -18,6 +18,11 @@ classdef BilinearForm < handle
%% methods
methods (Access=public)
mat = assemble(obj, fes);
+ mat = assemblePatchwise(obj, fes);
+ end
+
+ methods (Access=protected)
+ data = computeVolumeData(obj, fes);
end
methods (Static, Access=protected)
diff --git a/lib/assembly/@BilinearForm/assemble.m b/lib/assembly/@BilinearForm/assemble.m
index 400fc151e01c04a9186a0a8d9a08f8f3dacbaf48..7a33f904e97f751d8983841258ef6356cad3cf64 100644
--- a/lib/assembly/@BilinearForm/assemble.m
+++ b/lib/assembly/@BilinearForm/assemble.m
@@ -10,35 +10,8 @@ arguments
obj
fes FeSpace
end
-%% setup
-if all([isempty(obj.a), isempty(obj.b), isempty(obj.c), isempty(obj.robin)])
- error('BilinearForm:NoCoefficients', 'No coefficients are given.')
-end
-
-%% integrate element data
-mat = 0;
-phi = TestFunction(fes);
-Dphi = TestFunctionGradient(fes);
-
-% diffusion
-if ~isempty(obj.a)
- f = CompositeFunction(@BilinearForm.diffusionPart, obj.a, Dphi);
- mat = mat + integrateElement(f, obj.qra);
-end
-
-% convection
-if ~isempty(obj.b)
- f = CompositeFunction(@BilinearForm.convectionPart, obj.b, Dphi, phi);
- mat = mat + integrateElement(f, obj.qrb);
-end
-
-% reaction
-if ~isempty(obj.c)
- f = CompositeFunction(@BilinearForm.reactionPart, obj.c, phi);
- mat = mat + integrateElement(f, obj.qrc);
-end
-
%% construct sparse matrix
+mat = computeVolumeData(obj, fes);
dofs = getDofs(fes);
if ~isequal(mat, 0)
@@ -52,6 +25,7 @@ end
%% integrate and add boundary data
% Robin data
if ~isempty(obj.robin)
+ phi = TestFunction(fes);
f = CompositeFunction(@BilinearForm.robinPart, obj.robin, phi);
idx = getCombinedBndEdges(fes.mesh, fes.bnd.robin);
[I, J] = getLocalDofs(size(dofs.edge2Dofs, Dim.Vector));
diff --git a/lib/assembly/@BilinearForm/assemblePatchwise.m b/lib/assembly/@BilinearForm/assemblePatchwise.m
new file mode 100644
index 0000000000000000000000000000000000000000..036f370b2328dc311010a86b990403266725dea2
--- /dev/null
+++ b/lib/assembly/@BilinearForm/assemblePatchwise.m
@@ -0,0 +1,23 @@
+% assemblePatchwise Assembles matrix of bilinear form patchwise.
+%
+% A = assemblePatchwise(obj, fes) assembles the data of the given
+% BilinearForm for a FeSpace fes and returns a PatchwiseMatrix A.
+%
+% See also: BilinearForm
+
+function mat = assemblePatchwise(obj, fes, patches)
+arguments
+ obj
+ fes FeSpace
+ patches {mustBeIndexVector} = ':'
+end
+
+if ~isempty(obj.robin)
+ error('Patchwise assembly not implemented for Robin boundary!')
+end
+
+% in the long run, this should work just with non-assembled data:
+data = computeVolumeData(obj, fes);
+mat = PatchwiseMatrix(fes, data, patches);
+
+end
\ No newline at end of file
diff --git a/lib/assembly/@BilinearForm/computeVolumeData.m b/lib/assembly/@BilinearForm/computeVolumeData.m
new file mode 100644
index 0000000000000000000000000000000000000000..d13e55444301092bdfbaea550391ba11608ac4c5
--- /dev/null
+++ b/lib/assembly/@BilinearForm/computeVolumeData.m
@@ -0,0 +1,41 @@
+% computeVolumeData Computes raw data of bilinear form for volume terms.
+%
+% data = computeData(obj, fes) computes the elementwise data of the given
+% BilinearForm for a FeSpace fes for every non-trivial dof-pairing.
+%
+% See also: BilinearForm
+
+function data = computeVolumeData(obj, fes)
+arguments
+ obj
+ fes FeSpace
+end
+
+if all([isempty(obj.a), isempty(obj.b), isempty(obj.c), isempty(obj.robin)])
+ error('BilinearForm:NoCoefficients', 'No coefficients are given.')
+end
+
+%% integrate element data
+data = 0;
+phi = TestFunction(fes);
+Dphi = TestFunctionGradient(fes);
+
+% diffusion
+if ~isempty(obj.a)
+ f = CompositeFunction(@BilinearForm.diffusionPart, obj.a, Dphi);
+ data = data + integrateElement(f, obj.qra);
+end
+
+% convection
+if ~isempty(obj.b)
+ f = CompositeFunction(@BilinearForm.convectionPart, obj.b, Dphi, phi);
+ data = data + integrateElement(f, obj.qrb);
+end
+
+% reaction
+if ~isempty(obj.c)
+ f = CompositeFunction(@BilinearForm.reactionPart, obj.c, phi);
+ data = data + integrateElement(f, obj.qrc);
+end
+
+end
diff --git a/lib/assembly/@PatchwiseMatrix/PatchwiseMatrix.m b/lib/assembly/@PatchwiseMatrix/PatchwiseMatrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..8f68725a98dc187dcd0b481f9fd404f03491f519
--- /dev/null
+++ b/lib/assembly/@PatchwiseMatrix/PatchwiseMatrix.m
@@ -0,0 +1,63 @@
+% PatchwiseMatrix (handle class) Patchwise matrix of bilinear form.
+
+classdef PatchwiseMatrix < handle
+ properties (GetAccess=public, SetAccess=protected)
+ fes
+ activePatches
+ end
+
+ properties (Access=protected)
+ patchDofs
+ patchElements
+ patchesAsFreeDofs
+ patchwiseChol
+ end
+
+ properties (Dependent)
+ nPatches
+ end
+
+ methods (Access=public)
+ function obj = PatchwiseMatrix(fes, data, patches)
+ arguments
+ fes FeSpace
+ data
+ patches {mustBeIndexVector}
+ end
+
+ obj.fes = fes;
+ [obj.patchDofs, obj.patchElements] = assemblePatchDofs(obj);
+ idx = 1:obj.nPatches;
+ obj.activePatches = reshape(idx(patches), 1, []);
+ obj.buildLocalMatrixFactorization(data);
+ end
+
+ function dofs = getPatchDofs(obj, idx)
+ arguments
+ obj
+ idx {mustBeIndexVector}
+ end
+ dofs = vertcat(obj.patchDofs{idx});
+ end
+
+ function elements = getPatchElements(obj, idx)
+ arguments
+ obj
+ idx {mustBeIndexVector}
+ end
+ elements = vertcat(obj.patchElements{idx});
+ end
+ end
+
+ methods
+ function nPatches = get.nPatches(obj)
+ nPatches = numel(obj.patchDofs);
+ end
+ end
+
+ methods (Access=protected)
+ [patchDofs, patchElements] = assemblePatchDofs(obj)
+ buildLocalMatrixFactorization(obj, data)
+ end
+end
+
diff --git a/lib/assembly/@PatchwiseMatrix/assemblePatchDofs.m b/lib/assembly/@PatchwiseMatrix/assemblePatchDofs.m
new file mode 100644
index 0000000000000000000000000000000000000000..c4cd1a8f9145d8bcafd95c70f8bee17da3b3a4b5
--- /dev/null
+++ b/lib/assembly/@PatchwiseMatrix/assemblePatchDofs.m
@@ -0,0 +1,83 @@
+% assemblePatchDofs Assembles lists of unique inner dofs for each patch in
+% underlying mesh.
+%
+% assemblePatchDofs(obj)
+
+function [patchDofs, patchElements] = assemblePatchDofs(obj)
+%% preallocate arrays
+fes = obj.fes;
+mesh = fes.mesh;
+
+%% compute patch membership of (free) vertices/edges/elements
+dirichletEdges = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
+vertexNotOnDirichletBnd = ~ismember(1:mesh.nCoordinates, mesh.edges(:,dirichletEdges));
+
+free = computeFreeEdges(mesh, dirichletEdges);
+[patchEdges, nEdgesPerPatch] = computePatchMembership(mesh.edges(:,free));
+patchEdges = free(patchEdges);
+
+[patchElements, nElementsPerPatch] = computePatchMembership(mesh.elements);
+
+%% combine vertex/edge/element dofs for each patch
+vertexDofs = asVector(createVertexDofs(fes, vertexNotOnDirichletBnd));
+edgeDofs = asVector(createInnerEdgeDofs(fes, patchEdges));
+elementDofs = asVector(createInnerElementDofs(fes, patchElements));
+
+nLocalDofs = getDofConnectivity(fes.finiteElement);
+[patchDofs, nDofsPerPatch] = interleaveGroupedArrays(...
+ vertexDofs, nLocalDofs(1)*double(vertexNotOnDirichletBnd), ...
+ edgeDofs, nLocalDofs(2)*nEdgesPerPatch', ...
+ elementDofs, nLocalDofs(3)*nElementsPerPatch');
+
+patchDofs = mat2cell(patchDofs', nDofsPerPatch);
+patchElements = mat2cell(patchElements, nElementsPerPatch);
+
+end
+
+%% auxiliary functions
+
+% compute free edges without setdiff -> linear complexity
+function freeEdges = computeFreeEdges(mesh, dirichlet)
+ freeEdges = true(mesh.nEdges,1);
+ freeEdges(dirichlet) = false;
+ freeEdges = find(freeEdges)';
+end
+
+% requires one global sort of connectivity arrays -> log-linear complexity
+function [patches, count] = computePatchMembership(objects)
+ n = size(objects, 1);
+ N = size(objects, 2);
+
+ count = accumarray(objects(:), 1);
+ ptr = cumsum([1; count]);
+ patches = zeros(n*N, 1);
+
+ for i = 1:N
+ for k = 1:n
+ idx = objects(k,i);
+ patches(ptr(idx)) = i;
+ ptr(idx) = ptr(idx) + 1;
+ end
+ end
+end
+
+function [totalArray, totalGroupSize] = interleaveGroupedArrays(array, groupSize)
+arguments (Repeating)
+ array
+ groupSize
+end
+
+totalGroupPointer = cumsum([1; reshape(vertcat(groupSize{:}), [], 1)])';
+totalArray = zeros(1, totalGroupPointer(end)-1);
+totalGroupSize = sum(vertcat(groupSize{:}), 1);
+
+N = numel(groupSize);
+for n = 1:N
+ % add vector [0 1 2 0 1 0 1 2 3 ...] to group starting locations to
+ % make space for all group members
+ offset = totalGroupPointer(n:N:(end-1)) - cumsum([0, groupSize{n}(1:(end-1))]);
+ idx = (1:numel(array{n})) + repelem(offset, groupSize{n}) - 1;
+ totalArray(idx) = array{n};
+end
+
+end
diff --git a/lib/assembly/@PatchwiseMatrix/buildLocalMatrixFactorization.m b/lib/assembly/@PatchwiseMatrix/buildLocalMatrixFactorization.m
new file mode 100644
index 0000000000000000000000000000000000000000..e17b70cb22e4eec8b8eaf7f5edfd97ee29ecc8bd
--- /dev/null
+++ b/lib/assembly/@PatchwiseMatrix/buildLocalMatrixFactorization.m
@@ -0,0 +1,38 @@
+% buildLocalMatrixFactorization Patchwise assembly of matrix and cholesky
+% factorization.
+
+function buildLocalMatrixFactorization(obj, data)
+
+% since this operates only on free dofs, store global numbering
+freeDofs = getFreeDofs(obj.fes);
+dofs = getDofs(obj.fes);
+global2freeDofs = zeros(dofs.nDofs, 1);
+global2freeDofs(freeDofs) = 1:numel(freeDofs);
+obj.patchesAsFreeDofs = cellfun(@(x) global2freeDofs(x), obj.patchDofs, 'UniformOutput', false);
+nDofsPerElement = size(dofs.element2Dofs, 1);
+
+% store cholesky decomposition of local matrices associated to patches
+obj.patchwiseChol = cell(obj.nPatches, 1);
+for k = obj.activePatches
+ % get free dofs on patch and all dofs on patch
+ pDofs = obj.patchDofs{k};
+ n = numel(pDofs);
+ pElements = obj.patchElements{k};
+ allPatchDofs = dofs.element2Dofs(:,pElements);
+
+ % find local numbering of free dofs and relation between all and free dofs
+ [isNeeded, localNumbering] = ismember(allPatchDofs, pDofs);
+ isNeeded = reshape(isNeeded, size(isNeeded,1), 1, size(isNeeded,2));
+ idx = find(and(isNeeded, pagetranspose(isNeeded)));
+
+ % assemble arrays for sparse matrix with all dofs on patch ...
+ I = repmat(localNumbering, nDofsPerElement, 1);
+ J = repelem(localNumbering, nDofsPerElement, 1);
+ localData = data(:,pElements);
+
+ % ... and assemble/decompose only the parts that correspond to free dofs
+ obj.patchwiseChol{k} = ...
+ chol(accumarray([I(idx), J(idx)], localData(idx), [n, n]));
+end
+
+end
\ No newline at end of file
diff --git a/lib/assembly/@PatchwiseMatrix/mldivide.m b/lib/assembly/@PatchwiseMatrix/mldivide.m
new file mode 100644
index 0000000000000000000000000000000000000000..b8a157ec68a71e99a40422befff3d5bfdc444bc3
--- /dev/null
+++ b/lib/assembly/@PatchwiseMatrix/mldivide.m
@@ -0,0 +1,21 @@
+% mldivide Solve linear System A * x = y for patchwise matrix A. This means
+% solving local problems on patches associated to vertices
+% (depending on level) and add the patchwise solutions
+% together. mldivide(A, y) can also be called by A \ y.
+
+function x = mldivide(obj, y)
+
+x = zeros(size(y));
+lower = struct('UT', true, 'TRANSA', true);
+upper = struct('UT', true);
+
+for k = obj.activePatches
+ U = obj.patchwiseChol{k};
+ idx = obj.patchesAsFreeDofs{k};
+
+ % solve local patch problems (update is additive)
+ update = linsolve(U, linsolve(U, y(idx,:), lower), upper);
+ x(idx,:) = x(idx,:) + update;
+end
+
+end
diff --git a/lib/functions/@Evaluable/plot.m b/lib/functions/@Evaluable/plot.m
index 846bcdbad41b2d9775af7d91458d34b461b9acf6..995d7a9c25ab0d7ccea95673f4d6deb69fd52664 100644
--- a/lib/functions/@Evaluable/plot.m
+++ b/lib/functions/@Evaluable/plot.m
@@ -26,7 +26,10 @@ end
data = eval(obj, Barycentric2D([1;0;0]), 1);
nComponents = size(data, Dim.Vector);
ax = cell(nComponents, 1);
-for i = 1:nComponents
+ax{1} = gca;
+cla(ax{1})
+% TODO better handling of several components
+for i = 2:nComponents
ax{i} = axes(figure());
end
diff --git a/lib/mesh/@Mesh/Mesh.m b/lib/mesh/@Mesh/Mesh.m
index 48544b7906a04d13b775b2d01df769d33c1c660e..45fe85963d81fdba976e9c5ccac0a76d07898d22 100644
--- a/lib/mesh/@Mesh/Mesh.m
+++ b/lib/mesh/@Mesh/Mesh.m
@@ -100,6 +100,7 @@ classdef Mesh < handle
%% static methods
methods (Static)
obj = loadFromGeometry(geometryIdentifier)
+ obj = unitTriangle()
[edges, element2edges, flipEdges, boundary2edges] = computeEdgeInformation(elements, boundaries)
end
end
diff --git a/lib/mesh/@Mesh/unitTriangle.m b/lib/mesh/@Mesh/unitTriangle.m
new file mode 100644
index 0000000000000000000000000000000000000000..fbc0c23a0651786320034cf63e55d50699a2523e
--- /dev/null
+++ b/lib/mesh/@Mesh/unitTriangle.m
@@ -0,0 +1,19 @@
+% unitTriangle Construct Mesh object on unit triangle.
+%
+% mesh = Mesh.unitTriangle() constructs a mesh with one element and
+% homogeneous dirichlet boundary conditions on the unit triangle. In
+% particular, coordinates, elements, and edges are the same as in the
+% mesh constructed by mesh = Mesh.loadGeometry('unittriangle'). This
+% method, however, is faster.
+%
+% See also: Mesh
+
+function obj = unitTriangle()
+
+coordinates = [0, 1, 0; 0, 0, 1];
+elements = [2;3;1];
+boundaries = {[1, 2, 3; 2, 3, 1]};
+
+obj = Mesh(coordinates, elements, boundaries);
+
+end
\ No newline at end of file
diff --git a/lib/solvers/@IterativeSolver/IterativeSolver.m b/lib/solvers/@IterativeSolver/IterativeSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..be5b90b45bbca63d72b0cb500295acc8cb234787
--- /dev/null
+++ b/lib/solvers/@IterativeSolver/IterativeSolver.m
@@ -0,0 +1,133 @@
+% IterativeSolver (abstract handle class) Interface for iterative solvers.
+%
+% solver.setupSystemMatrix(A) sets up system Matrix for the linear system
+% A*x=b. The matrix A must be quadratic.
+%
+% solver.setupRhs(b, [x0]) sets up right-hand side for the linear system A*x=b
+% and initial guess x0 (0 per default). The right-hand side b can have
+% multiple columns, but b and x0 must have the same size.
+%
+% solver.solve() performs an automatic loop of solver.step() until
+% convergence in each column of the right-hand side is reached, checking
+% the stopping criterion after every step.
+%
+% solver.step() performs one step of the iterative solver on active
+% components of the right-hand side and automatoically keeps track of
+% iteration numbers.
+%
+% tf = solver.applyStoppingCriterion() checks stopping criterion and set
+% converged components to inactive.
+%
+% Abstract methods to implement:
+%
+% solver.computeUpdate() computes the update of active components of the
+% solution.
+%
+% tf = solver.isConverged(solver) returns state of convergence of the solver
+% for each column of the right-hand side.
+%
+% [solver, P] = IterativeSolver.choose(fes, blf, class, variant) convenience
+% factory method. Returns instances of IterativeSolver and Prolongation
+% accoding to the given parameters.
+
+
+classdef IterativeSolver < handle
+ %% properties
+ properties (Access=public)
+ maxIter (1,:) double {mustBePositive}
+ tol (1,:) double {mustBePositive}
+ end
+
+ properties (SetAccess=protected, GetAccess=public)
+ A
+ b
+ x
+ iterationCount
+ activeComponents
+ end
+
+ methods (Abstract, Access=public)
+ isConverged(obj)
+ end
+
+ methods (Abstract, Access=protected)
+ computeUpdate(obj)
+ end
+
+ %% generic routines for iterative solvers
+ % overwrite to include additional steps
+ methods (Access=public)
+ function obj = IterativeSolver()
+ obj.maxIter = 100;
+ obj.tol = 1e-8;
+ end
+
+ function setupSystemMatrix(obj, A)
+ arguments
+ obj
+ A (:,:) double {IterativeSolver.mustBeQuadratic}
+ end
+ obj.A = A;
+ end
+
+ function setupRhs(obj, b, x0)
+ arguments
+ obj
+ b (:,:) double
+ x0 (:,:) double {IterativeSolver.mustBeCompatible(b, x0)} = zeros(size(b));
+ end
+ assert(isequal(size(obj.A,2), size(b,1)), ...
+ 'Right-hand size and matrix must be of compatible size.')
+ obj.b = b;
+ obj.x = x0;
+ obj.iterationCount = zeros(1, size(b,2));
+ obj.activeComponents = true(1, size(b,2));
+ end
+
+ function x = solve(obj)
+ while ~all(applyStoppingCriterion(obj))
+ obj.step()
+ end
+
+ if any(obj.iterationCount >= obj.maxIter)
+ warning('Maximal number of iterations exhausted! maxIter = %d', obj.maxIter)
+ end
+
+ x = obj.x;
+ end
+
+ function step(obj)
+ obj.computeUpdate();
+ obj.iterationCount(obj.activeComponents) = obj.iterationCount(obj.activeComponents) + 1;
+ end
+
+ function tf = applyStoppingCriterion(obj)
+ tf = isConverged(obj);
+ obj.activeComponents = obj.activeComponents & ~tf;
+ end
+ end
+
+ %% convenience factory function
+ methods (Static, Access=public)
+ [solver, P] = choose(fes, blf, class, variant)
+ end
+
+ %% validation functions to use within this class
+ methods (Static, Access=protected)
+ function mustBeQuadratic(A)
+ if ~(size(A,1) == size(A,2))
+ eidType = 'Solver:matrixNotQuadratic';
+ msgType = 'Matrix must be quadratic.';
+ throwAsCaller(MException(eidType,msgType))
+ end
+ end
+
+ function mustBeCompatible(b, x0)
+ if ~isempty(x0) && ~isequal(size(b), size(x0))
+ eidType = 'Solver:dataNotCompatible';
+ msgType = 'Size of RHS and initial guess must be compatible.';
+ throwAsCaller(MException(eidType,msgType))
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/@IterativeSolver/choose.m b/lib/solvers/@IterativeSolver/choose.m
new file mode 100644
index 0000000000000000000000000000000000000000..ef5d4fa6dbc32592f5c496f354ff839080dbd3a7
--- /dev/null
+++ b/lib/solvers/@IterativeSolver/choose.m
@@ -0,0 +1,91 @@
+% choose Return suitable solver (instance of subclass of IterativeSolver) and
+% suitable Prolongation P.
+%
+% solver = choose(fes, blf, class) returns a solver of class 'class' for given
+% finite element space fes and bilinear form blf. For each class, sensible
+% default variants are chosen.
+%
+% solver = choose(fes, blf, class, variant) returns a solver of class 'class'
+% of specific variant 'variant'.
+%
+% [solver, P] = choose(__) additionally to the solver, returns a suitable
+% Prolongation P, that can be used to prolongate the solution of the
+% solver to the next finer mesh.
+
+function [solver, P] = choose(fes, blf, class, variant)
+arguments
+ fes FeSpace
+ blf BilinearForm
+ class (1,1) string {mustBeMember(class, ["cg", "pcg", "multigrid", "direct"])}
+ variant (1,1) string {mustBeMember(variant, [ "", ...
+ "iChol", "jacobi", ...
+ "additiveSchwarzLowOrder", "additiveSchwarzHighOrder" ...
+ "lowOrderVcycle", "highOrderVcycle"])} = ""
+end
+
+order = fes.finiteElement.order;
+P = Prolongation.chooseFor(fes);
+
+switch class
+ % non-preconditioned CG
+ case "cg"
+ solver = PcgSolver(NoPreconditioner());
+
+ % preconditioned CG family
+ case "pcg"
+ switch variant
+ case "iChol"
+ preconditioner = IncompleteCholesky();
+ case "jacobi"
+ if order == 1
+ preconditioner = DiagonalJacobi();
+ else
+ preconditioner = BlockJacobi(fes, blf);
+ end
+ case {"", "additiveSchwarzLowOrder"}
+ if order == 1
+ preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
+ else
+ preconditioner = OptimalMLAdditiveSchwarz(JacobiLowOrderCascade(fes, blf));
+ end
+ case "additiveSchwarzHighOrder"
+ if order == 1
+ preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
+ else
+ preconditioner = OptimalMLAdditiveSchwarz(JacobiHighOrderCascade(fes, blf, P));
+ end
+ otherwise
+ error('No PCG variant %s!', variant)
+ end
+ solver = PcgSolver(preconditioner);
+
+ % geometric , Pmultigrid family
+ case "multigrid"
+ switch variant
+ case {"", "lowOrderVcycle"}
+ if order == 1
+ smoother = P1JacobiCascade(fes, blf, P);
+ else
+ smoother = JacobiLowOrderCascade(fes, blf);
+ end
+ case "highOrderVcycle"
+ if order == 1
+ smoother = P1JacobiCascade(fes, blf, P);
+ else
+ smoother = JacobiHighOrderCascade(fes, blf, P);
+ end
+ otherwise
+ error('No multigrid variant %s!', variant)
+ end
+ solver = OptimalVcycleMultigridSolver(smoother);
+
+ % direct solve (for testing purposes)
+ case "direct"
+ solver = DirectSolver();
+
+ % default
+ otherwise
+ error('No iterative solver of class %s!', class)
+end
+
+end
diff --git a/lib/solvers/CgSolver.m b/lib/solvers/CgSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..acfb831c4c5c86c9828e73605b5149eb0b51af9d
--- /dev/null
+++ b/lib/solvers/CgSolver.m
@@ -0,0 +1,12 @@
+% CgSolver (subclass of PcgSolver) Solves linear equations iteratively using
+% the CG method without preconditioner.
+
+classdef CgSolver < PcgSolver
+ %% methods
+ methods (Access=public)
+ % preconditioner = identity
+ function Cx = preconditionAction(obj, x)
+ Cx = x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/DirectSolver.m b/lib/solvers/DirectSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..91ca5446ade055b0a071bc3843a1fddee82e9e28
--- /dev/null
+++ b/lib/solvers/DirectSolver.m
@@ -0,0 +1,21 @@
+% DirectSolver (subclass of IterativeSolver) Solves linear equations
+% directly. This is mainly intended for testing purposes!
+%
+% See also: IterativeSolver
+
+classdef DirectSolver < IterativeSolver
+ %% implement abstract superclass methods
+ methods (Access=public)
+ function tf = isConverged(obj)
+ % converges in the first step
+ tf = (obj.iterationCount > 0);
+ end
+ end
+
+ methods (Access=protected)
+ function computeUpdate(obj)
+ % direct solve
+ obj.x = obj.A \ obj.b;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/OptimalVcycleMultigridSolver.m b/lib/solvers/OptimalVcycleMultigridSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..496c2a73082af19af51919bf7ebb1c0ae04f8572
--- /dev/null
+++ b/lib/solvers/OptimalVcycleMultigridSolver.m
@@ -0,0 +1,144 @@
+% OptimalVcycleMultigridSolver (subclass of IterativeSolver) Solves linear
+% equations iteratively, using a geometric multigrid solver with specified
+% MultiLevelSmoother.
+% One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
+%
+% To obtain the solver described in [Innerberger, Miraçi, Praetorius,
+% Streitberger; Optimal computational costs of AFEM with optimal local
+% hp-robust multigrid solver], choose the smoother as DiagonalJacobi (p=1),
+% LowOrderBlockJacobi (p > 1, 1->1->p iteration), or HighOrderBlockJacobi
+% (p > 1, 1->p->p iteration).
+% This solver also comes with built-in estimator for the algebraic error,
+% which is also returned by this multigrid solver.
+%
+% See also IterativeSolver, MultiLevelSmoother, DiagonalJacobi,
+% LowOrderBlockJacobi, HighOrderBlockJacobi
+
+classdef OptimalVcycleMultigridSolver < IterativeSolver
+ %% properties
+ properties (SetAccess=protected, GetAccess=public)
+ residualCNorm
+ algEstimator
+ residual
+ Cresidual
+ end
+
+ properties (Dependent, Access=public)
+ nLevels
+ end
+
+ properties (Access=protected)
+ smoother
+ projectedMatrix
+ lhsNorm
+ end
+
+ methods
+ function obj = OptimalVcycleMultigridSolver(smoother)
+ obj = obj@IterativeSolver();
+ obj.smoother = smoother;
+ end
+
+ function n = get.nLevels(obj)
+ n = obj.smoother.nLevels;
+ end
+
+ function setupSystemMatrix(obj, A, varargin)
+ % override superclass method, because ALL matrices are needed
+ setupSystemMatrix@IterativeSolver(obj, A, varargin{:});
+
+ projectedA = obj.smoother.setup(A, varargin{:});
+ obj.projectedMatrix{obj.nLevels} = projectedA;
+ end
+
+ function setupRhs(obj, b, varargin)
+ setupRhs@IterativeSolver(obj, b, varargin{:});
+ % initialize residual & hierarchy
+ obj.residual = b - obj.A * obj.x;
+ [obj.Cresidual, obj.algEstimator] = obj.Vcycle(obj.residual);
+ obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, Dim.Vector));
+ obj.lhsNorm = vecnorm(b, 2, Dim.Vector);
+ end
+
+ % stopping criterion
+ function tf = isConverged(obj)
+ tf = ((obj.iterationCount >= obj.maxIter) ...
+ | (obj.residualCNorm < obj.tol) ...
+ | (obj.residualCNorm ./ obj.lhsNorm < obj.tol));
+ end
+ end
+
+ methods (Access=protected)
+ % one step
+ function computeUpdate(obj)
+ % only iterate on active components
+ idx = find(obj.activeComponents);
+
+ % update solution & residual
+ obj.x(:,idx) = obj.x(:,idx) + obj.Cresidual(:,idx);
+ obj.residual(:,idx) = obj.residual(:,idx) - obj.A * obj.Cresidual(:,idx);
+
+ % residual & update error correction
+ [obj.Cresidual(:,idx), obj.algEstimator(idx)] = obj.Vcycle(obj.residual(:,idx));
+ obj.residualCNorm(:,idx) = sqrt(dot(obj.residual(:,idx), obj.Cresidual(:,idx), Dim.Vector));
+ end
+
+ function [Cx, algEstimator] = Vcycle(obj, x)
+ assert(isequal(size(x, 1), size(obj.A, 1)), ...
+ 'Setup for multilevel iteration not complete!')
+
+ % if there is only one coarse level: exact solve
+ L = obj.nLevels;
+ if L == 1
+ Cx = obj.A \ x;
+ algEstimator = sqrt(dot(Cx, obj.A * Cx, Dim.Vector));
+ return
+ end
+
+ % descending cascade: no smoothing
+ res{L} = x;
+ for k = L:-1:2
+ res{k-1} = obj.smoother.restrict(res{k}, k);
+ end
+
+ % exact solve on coarsest level to compute accumulative lifting of x
+ Cx = obj.projectedMatrix{1} \ res{1};
+ eta2 = dot(Cx, obj.projectedMatrix{1} * Cx, Dim.Vector);
+
+ % ascending cascade: (local) smoothing and compute error estimator
+ for k = 2:(L-1)
+ Cx = obj.smoother.prolongate(Cx, k);
+ updatedResidual = res{k} - obj.projectedMatrix{k} * Cx;
+ rho = obj.smoother.smooth(updatedResidual, k);
+ [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.projectedMatrix{k}, updatedResidual, rho);
+ Cx = Cx + xUpdate;
+ eta2 = eta2 + etaUpdate2;
+ end
+
+ % final smoothing and residual update (with non-projected matrix)
+ Cx = obj.smoother.prolongate(Cx, L);
+ updatedResidual = res{L} - obj.A * Cx;
+ rho = obj.smoother.smooth(updatedResidual, L);
+ [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.A, updatedResidual, rho);
+ Cx = Cx + xUpdate;
+ algEstimator = sqrt(eta2 + etaUpdate2);
+ end
+ end
+end
+
+%% auxiliary functions
+% error correction with optimal stepsize
+function [etaUpdate2, resUpdate] = computeOptimalUpdate(A, res, rho)
+ rhoArho = dot(rho, A*rho, Dim.Vector);
+ if max(abs(rho)) < eps
+ lambda = 1;
+ else
+ lambda = dot(res, rho, Dim.Vector) ./ rhoArho;
+ end
+ resUpdate = lambda .* rho;
+ etaUpdate2 = rhoArho .* lambda.^2;
+
+ if any(lambda > 3)
+ warning('MG step-sizes no longer bound by d+1. Optimality of step size cannot be guaranteed!')
+ end
+end
diff --git a/lib/solvers/PcgSolver.m b/lib/solvers/PcgSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..1a278c82dda3a969218f6e1d6f730baa11c4bdcc
--- /dev/null
+++ b/lib/solvers/PcgSolver.m
@@ -0,0 +1,86 @@
+% PcgSolver (abstract subclass of IterativeSolver) Solves linear equations
+% iteratively by a preconditioned conjugate gradient (PCG) method.
+%
+% See also: IterativeSolver
+
+classdef PcgSolver < IterativeSolver
+ %% properties
+ properties (SetAccess=protected, GetAccess=public)
+ residualCNorm
+ end
+
+ properties (Access=protected)
+ C
+ residual
+ Cresidual
+ searchDirection
+ normb
+ end
+
+ %% constructor
+ methods (Access=public)
+ function obj = PcgSolver(preconditioner)
+ arguments
+ preconditioner (1,1) Preconditioner
+ end
+
+ obj = obj@IterativeSolver();
+ obj.C = preconditioner;
+ end
+ end
+
+ %% extend superclass methods
+ methods (Access=public)
+ function setupSystemMatrix(obj, A)
+ setupSystemMatrix@IterativeSolver(obj, A);
+ obj.C.setup(A);
+ end
+
+ function setupRhs(obj, b, varargin)
+ setupRhs@IterativeSolver(obj, b, varargin{:});
+
+ % initialize residual & search direction
+ obj.residual = b - obj.A * obj.x;
+ obj.Cresidual = obj.C.apply(obj.residual);
+ obj.searchDirection = obj.Cresidual;
+ obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, 1));
+ obj.normb = vecnorm(b, 2, 1);
+ end
+ end
+
+ %% implement abstract superclass methods
+ methods (Access=public)
+ function tf = isConverged(obj)
+ tf = ((obj.iterationCount >= obj.maxIter) ...
+ | (obj.residualCNorm < obj.tol) ...
+ | (obj.residualCNorm ./ obj.normb < obj.tol));
+ end
+ end
+
+ methods (Access=protected)
+ function computeUpdate(obj)
+ % only iterate on active components
+ idx = find(obj.activeComponents);
+
+ % update solution
+ AsearchDirection = obj.A * obj.searchDirection(:,idx);
+ dAd = dot(obj.searchDirection(:,idx), AsearchDirection, 1);
+ if dAd < eps
+ alpha = 1;
+ else
+ alpha = obj.residualCNorm(:,idx).^2 ./ dAd;
+ end
+ obj.x(:,idx) = obj.x(:,idx) + alpha .* obj.searchDirection(:,idx);
+
+ % update residual
+ obj.residual(:,idx) = obj.residual(:,idx) - alpha .* AsearchDirection;
+ obj.Cresidual(:,idx) = obj.C.apply(obj.residual(:,idx));
+ residualCNormOld = obj.residualCNorm(:,idx);
+ obj.residualCNorm(:,idx) = sqrt(dot(obj.residual(:,idx), obj.Cresidual(:,idx), 1));
+
+ % update search direction
+ beta = (obj.residualCNorm(:,idx) ./ residualCNormOld).^2;
+ obj.searchDirection(:,idx) = obj.Cresidual(:,idx) + beta .* obj.searchDirection(:,idx);
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/BlockJacobi.m b/lib/solvers/preconditioners/BlockJacobi.m
new file mode 100644
index 0000000000000000000000000000000000000000..d4bf3ec1cffaab76bcea9d48ed8d4f036e92a34c
--- /dev/null
+++ b/lib/solvers/preconditioners/BlockJacobi.m
@@ -0,0 +1,34 @@
+% BlockJacobi (subclass of Preconditioner) block-diagonal preconditioner.
+%
+% See also: Preconditioner, PcgSolver
+
+classdef BlockJacobi < Preconditioner
+ %% properties
+ properties (Access=private)
+ C
+ fes
+ blf
+ end
+
+ %% methods
+ methods (Access=public)
+ % preconditioner action: patchwise inverse of diagonal
+ function obj = BlockJacobi(fes, blf)
+ arguments
+ fes FeSpace
+ blf BilinearForm
+ end
+
+ obj.fes = fes;
+ obj.blf = blf;
+ end
+
+ function setup(obj, A)
+ obj.C = assemblePatchwise(obj.blf, obj.fes, ':');
+ end
+
+ function Cx = apply(obj, x)
+ Cx = obj.C \ x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/DiagonalJacobi.m b/lib/solvers/preconditioners/DiagonalJacobi.m
new file mode 100644
index 0000000000000000000000000000000000000000..9ade2832478bd8c204fcb6c3c99b89e5c1132500
--- /dev/null
+++ b/lib/solvers/preconditioners/DiagonalJacobi.m
@@ -0,0 +1,22 @@
+% DiagonalJacobi (subclass of Preconditioner) diagonal preconditioner.
+%
+% See also: Preconditioner, PcgSolver
+
+classdef DiagonalJacobi < Preconditioner
+ %% properties
+ properties (Access=private)
+ C
+ end
+
+ %% methods
+ methods (Access=public)
+ % preconditioner action: inverse of diagonal
+ function setup(obj, A)
+ obj.C = full(diag(A)).^(-1);
+ end
+
+ function Cx = apply(obj, x)
+ Cx = obj.C .* x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/IncompleteCholesky.m b/lib/solvers/preconditioners/IncompleteCholesky.m
new file mode 100644
index 0000000000000000000000000000000000000000..6a04a2577a04519a146d399cbe4ecbeb9b242583
--- /dev/null
+++ b/lib/solvers/preconditioners/IncompleteCholesky.m
@@ -0,0 +1,22 @@
+% IncompleteCholesky (subclass of Preconditioner) incomplete Cholesky
+% decomposition preconditioner.
+%
+% See also: Preconditioner, PcgSolver
+
+classdef IncompleteCholesky < Preconditioner
+ %% properties
+ properties (Access=private)
+ C
+ end
+
+ %% methods
+ methods (Access=public)
+ function setup(obj, A)
+ obj.C = ichol(A);
+ end
+
+ function Cx = apply(obj, x)
+ Cx = obj.C' \ (obj.C \ x);
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/NoPreconditioner.m b/lib/solvers/preconditioners/NoPreconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..fa413b5e0f21d730a0c570740656d09f603c7d5d
--- /dev/null
+++ b/lib/solvers/preconditioners/NoPreconditioner.m
@@ -0,0 +1,12 @@
+% NoPreconditioner (subclass of Preconditioner) identity preconditioner (=no-op).
+%
+% See also: Preconditioner, PcgSolver
+
+
+classdef NoPreconditioner < Preconditioner
+ methods (Access=public)
+ function Cx = apply(~, x)
+ Cx = x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m b/lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m
new file mode 100644
index 0000000000000000000000000000000000000000..c7499bdc6d56acac17c2531fd484ba6ac3c603b9
--- /dev/null
+++ b/lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m
@@ -0,0 +1,76 @@
+% OptimalMLAdditiveSchwarz (subclass of Preconditioner) optimal multilevel
+% additive Schwarz preconditioner with given smoother.
+%
+% See also: Preconditioner, PcgSolver, MultilevelSmoother
+
+classdef OptimalMLAdditiveSchwarz < Preconditioner
+ %% properties
+ properties (Access=protected)
+ smoother
+ Acoarse
+ projectedAcoarse
+ end
+
+ properties (Dependent)
+ nLevels
+ end
+
+ properties (Access=private)
+ listenerHandle
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = OptimalMLAdditiveSchwarz(smoother)
+ arguments
+ smoother MultilevelSmoother
+ end
+
+ obj = obj@Preconditioner();
+ obj.smoother = smoother;
+ end
+
+ function setup(obj, A)
+ if obj.nLevels == 0
+ obj.Acoarse = A;
+ obj.projectedAcoarse = obj.smoother.setup(A);
+ else
+ obj.smoother.setup(A);
+ end
+ end
+
+ function Cx = apply(obj, x)
+ assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
+
+ L = obj.nLevels;
+ rho = cell(L, 1);
+
+ if L == 1
+ Cx = obj.Acoarse \ x;
+ return
+ end
+
+ % descending cascade
+ residual = x;
+ for k = L:-1:2
+ rho{k} = obj.smoother.smooth(residual, k);
+ residual = obj.smoother.restrict(residual, k);
+ end
+
+ % exact solve on coarsest level
+ Cx = obj.projectedAcoarse \ residual;
+
+ % ascending cascade
+ for k = 2:L
+ Cx = obj.smoother.prolongate(Cx, k);
+ Cx = Cx + rho{k};
+ end
+ end
+ end
+
+ methods
+ function n = get.nLevels(obj)
+ n = obj.smoother.nLevels;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/preconditioners/Preconditioner.m b/lib/solvers/preconditioners/Preconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..88833d8618a7f99c5c25dc8e83edf8e46f840a16
--- /dev/null
+++ b/lib/solvers/preconditioners/Preconditioner.m
@@ -0,0 +1,23 @@
+% Preconditioner (abstract handle class) Interface for preconditioners.
+%
+% pc.setup(A, ...) hook to set up the preconditioner for the matrix A.
+% Default implementation does nothing.
+%
+% Abstract methods to implement:
+%
+% Cx = pc.apply(x) applies the preconditioner to the vector x.
+%
+% See also: PcgSolver
+
+
+classdef Preconditioner < handle
+ methods (Access=public)
+ function setup(~, ~)
+ % Default implementation does nothing.
+ end
+ end
+
+ methods (Abstract, Access=public)
+ apply(obj, x)
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/JacobiHighOrderCascade.m b/lib/solvers/smoothers/JacobiHighOrderCascade.m
new file mode 100644
index 0000000000000000000000000000000000000000..d58ab6995e87d15155a6a8aac5e58fe7a68779a2
--- /dev/null
+++ b/lib/solvers/smoothers/JacobiHighOrderCascade.m
@@ -0,0 +1,85 @@
+% JacobiHighOrderCascade (subclass of MultilevelSmoother) multilevel
+% Jacobi smoother for higher order finite elements: smooth locally with
+% patchwise stiffness matrix on changed patches on every level, except for
+% coarsest level, where global p1-solving is done. Because of the lifting from
+% p1 to higher order between the coarsest and the second level, all patches
+% are assumed to be changed on the second level
+%
+% Reference: [Innerberger, Miraçi, Praetorius, Streitberger; Optimal
+% computational costs of AFEM with optimal local hp-robust multigrid solver].
+%
+% See also: MultilevelSmoother, OptimalVcycleMultigridSolver
+
+
+classdef JacobiHighOrderCascade < MultilevelSmoother
+ properties (Access=protected)
+ P
+ loFes
+ patchwiseA
+ inclusionMatrix
+ intergridMatrix
+ freeDofs
+ freeDofsOld
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = JacobiHighOrderCascade(fes, blf, P)
+ obj = obj@MultilevelSmoother(fes, blf);
+
+ assert(fes.finiteElement.order > 1, ...
+ 'JacobiHighOrderCascadeSmoother only works for higher order finite elements.')
+
+ obj.P = P;
+ obj.loFes = FeSpace(fes.mesh, LowestOrderH1Fe(), 'dirichlet', fes.bnd.dirichlet);
+ end
+
+ function nonInvertedSmootherMatrix = setup(obj, A)
+ setup@MultilevelSmoother(obj, A);
+
+ L = obj.nLevels;
+ obj.freeDofsOld = obj.freeDofs;
+ obj.freeDofs = getFreeDofs(obj.fes);
+ freeVertices = getFreeDofs(obj.loFes);
+
+ hoFes = obj.fes;
+ if L == 1
+ % set up FE inclusion and lowest order matrix on coarsest level
+ obj.inclusionMatrix = SpaceProlongation(obj.loFes, hoFes) ...
+ .matrix(getFreeDofs(obj.loFes), getFreeDofs(hoFes));
+ nonInvertedSmootherMatrix = assemble(obj.blf, obj.loFes);
+ nonInvertedSmootherMatrix = nonInvertedSmootherMatrix(freeVertices, freeVertices);
+ else
+ % note: patchwiseA is not already indexed by free vertices -> use freeVertices here
+ obj.changedPatches{L} = freeVertices(obj.changedPatches{L}(freeVertices));
+ obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
+ nonInvertedSmootherMatrix = A;
+ if L == 2
+ obj.patchwiseA{L} = assemblePatchwise(obj.blf, hoFes, ':');
+ else
+ obj.patchwiseA{L} = assemblePatchwise(obj.blf, hoFes, obj.changedPatches{L});
+ end
+ end
+ end
+
+ function Cx = smooth(obj, x, k)
+ % higher order global (k=2) or local (else) patchwise smooting
+ Cx = obj.patchwiseA{k} \ x;
+ end
+
+ % order of FESpace per level: 1 -> p -> ... -> p
+ function Px = prolongate(obj, x, k)
+ if k == 2
+ x = obj.inclusionMatrix' * x;
+ end
+ Px = obj.intergridMatrix{k} * x;
+ end
+
+ function Px = restrict(obj, x, k)
+ Px = obj.intergridMatrix{k}' * x;
+ if k == 2
+ Px = obj.inclusionMatrix * Px;
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/JacobiLowOrderCascade.m b/lib/solvers/smoothers/JacobiLowOrderCascade.m
new file mode 100644
index 0000000000000000000000000000000000000000..b0b9e5545ee9e79ae3e8420dedf94b822bca81f3
--- /dev/null
+++ b/lib/solvers/smoothers/JacobiLowOrderCascade.m
@@ -0,0 +1,90 @@
+% JacobiLowOrderCascade (subclass of MultilevelSmoother) multilevel
+% Jacobi smoother for higher order finite elements: smooth locally with
+% diagonal of P1 stiffness matrix on changed patches on every level, except
+% for finest level, where global patchwise higher order smoothing is done.
+%
+% Reference: [Innerberger, Miraçi, Praetorius, Streitberger; Optimal
+% computational costs of AFEM with optimal local hp-robust multigrid solver].
+%
+% See also: MultilevelSmoother, OptimalVcycleMultigridSolver
+
+
+classdef JacobiLowOrderCascade < MultilevelSmoother
+ properties (Access=protected)
+ P
+ loFes
+ p1Smoother
+ hoSmoother
+ inclusionMatrix
+ intergridMatrix
+ freeVertices
+ freeVerticesOld
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = JacobiLowOrderCascade(fes, blf)
+ obj = obj@MultilevelSmoother(fes, blf);
+
+ assert(fes.finiteElement.order > 1, ...
+ 'JacobiLowOrderCascadeSmoother only works for higher order finite elements.')
+
+ mesh = fes.mesh;
+ obj.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', fes.bnd.dirichlet);
+ obj.P = Prolongation.chooseFor(obj.loFes);
+
+ end
+
+ function nonInvertedSmootherMatrix = setup(obj, A)
+ setup@MultilevelSmoother(obj, A);
+
+ L = obj.nLevels;
+ obj.freeVerticesOld = obj.freeVertices;
+ obj.freeVertices = getFreeDofs(obj.loFes);
+
+ p1Matrix = assemble(obj.blf, obj.loFes);
+ p1Matrix = p1Matrix(obj.freeVertices, obj.freeVertices);
+ obj.p1Smoother{L} = full(diag(p1Matrix)).^(-1);
+ nonInvertedSmootherMatrix = p1Matrix;
+
+ if L >= 2
+ % note: p1 smoother is already indexed by free vertices -> use find here
+ obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
+ % intermediate level smoothers and transfer operators
+ obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
+ obj.p1Smoother{L} = obj.p1Smoother{L}(obj.changedPatches{L});
+ % finest level smoother and transfer operator
+ obj.hoSmoother = assemblePatchwise(obj.blf, obj.fes, ':');
+ obj.inclusionMatrix = SpaceProlongation(obj.loFes, obj.fes) ...
+ .matrix(getFreeDofs(obj.loFes), getFreeDofs(obj.fes));
+ end
+ end
+
+ function Cx = smooth(obj, x, k)
+ if k == obj.nLevels
+ % higher order global patchwise smooting
+ Cx = obj.hoSmoother \ x;
+ else
+ % local p1 patchwise smoothing
+ idx = obj.changedPatches{k};
+ Cx = zeros(size(x));
+ Cx(idx,:) = obj.p1Smoother{k} .* x(idx,:);
+ end
+ end
+
+ % order of FESpace per level: 1 -> ... -> 1 -> p
+ function Px = prolongate(obj, x, k)
+ Px = obj.intergridMatrix{k} * x;
+ if k == obj.nLevels
+ Px = obj.inclusionMatrix' * Px;
+ end
+ end
+
+ function Px = restrict(obj, x, k)
+ if k == obj.nLevels
+ x = obj.inclusionMatrix * x;
+ end
+ Px = obj.intergridMatrix{k}' * x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/MultilevelSmoother.m b/lib/solvers/smoothers/MultilevelSmoother.m
new file mode 100644
index 0000000000000000000000000000000000000000..ca7271db059c51bb0e23f923b20bd0e54110dc7c
--- /dev/null
+++ b/lib/solvers/smoothers/MultilevelSmoother.m
@@ -0,0 +1,84 @@
+% MultilevelSmoother (abstract handle class) Interface for smoothers to use in a
+% multigrid solver.
+%
+% Abstract methods to implement:
+%
+% smoother.setup(A) hook to set up the smoother for the matrix A.
+%
+% nonInvertedSmootherMatrix = smoother.setup(A) set up the smoother and also
+% return the non-inverted matrix, based on which the smoother is computed.
+%
+% Cx = smoother.smooth(x, k) applies the smoother to the vector x on level k.
+%
+% Px = smoother.prolongate(x, k) prolongate vector x from level k to level k+1.
+%
+% Px = smoother.restrict(x, k) restrict vector x from level k+1 to level k.
+%
+% See also: OptimalVcycleMultigridSolver
+
+
+classdef MultilevelSmoother < handle
+ %% properties
+ properties (GetAccess=public, SetAccess=protected)
+ nLevels
+ end
+
+ properties (Access=protected)
+ fes
+ blf
+ changedPatches
+ end
+
+ properties (Access=private)
+ listenerHandle
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = MultilevelSmoother(fes, blf)
+ arguments
+ fes FeSpace
+ blf BilinearForm
+ end
+
+ assert(isempty(blf.b), ...
+ 'Multilevel smoothers only tested for symmetric problems.')
+
+ obj.nLevels = 0;
+ obj.fes = fes;
+ obj.blf = blf;
+
+ mesh = fes.mesh;
+ obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
+ end
+
+ function nonInvertedSmootherMatrix = setup(obj, ~)
+ obj.nLevels = obj.nLevels + 1;
+
+ % should be implemented by subclasses
+ nonInvertedSmootherMatrix = [];
+ end
+ end
+
+ methods (Access=protected)
+ % callback function to be executed before mesh refinement:
+ % -> patches change for all vertices adjacent to refined edges and
+ % all new vertices
+ function getChangedPatches(obj, mesh, bisecData)
+ nCOld = mesh.nCoordinates;
+ nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
+ bisecData.nRefinedElements' * bisecData.nInnerNodes;
+ bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
+ obj.changedPatches{obj.nLevels+1} = false(nCNew, 1);
+ idx = [unique(bisectedEdgeNodes); ((nCOld+1):nCNew)'];
+ obj.changedPatches{obj.nLevels+1}(idx) = true;
+ end
+ end
+
+ %% abstract methods
+ methods (Abstract, Access=public)
+ smooth(obj, x, k)
+ prolongate(obj, x, k)
+ restrict(obj, x, k)
+ end
+end
\ No newline at end of file
diff --git a/lib/solvers/smoothers/P1JacobiCascade.m b/lib/solvers/smoothers/P1JacobiCascade.m
new file mode 100644
index 0000000000000000000000000000000000000000..93d15999bc92e42a7ebcb05955c717d93cdd4844
--- /dev/null
+++ b/lib/solvers/smoothers/P1JacobiCascade.m
@@ -0,0 +1,69 @@
+% P1JacobiCascade (subclass of MultilevelSmoother) multilevel Jacobi smoother
+% for lowest order finite elements: smooth with diagonal of stiffness matrix
+% on changed patches on every level.
+%
+% See also: MultilevelSmoother, OptimalVcycleMultigridSolver
+
+
+classdef P1JacobiCascade < MultilevelSmoother
+ properties (Access=protected)
+ P
+ inverseDiagonal
+ intergridMatrix
+ freeVertices
+ freeVerticesOld
+ end
+
+ %% event data
+ properties (Access=private)
+ listenerHandle
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = P1JacobiCascade(fes, blf, P)
+ arguments
+ fes
+ blf
+ P Prolongation
+ end
+
+ assert(fes.finiteElement.order == 1, ...
+ 'P1JacobiSmoother only works for lowest order finite elements.')
+
+ obj = obj@MultilevelSmoother(fes, blf);
+ obj.P = P;
+ end
+
+ function nonInvertedSmootherMatrix = setup(obj, A)
+ setup@MultilevelSmoother(obj, A);
+
+ L = obj.nLevels;
+ obj.freeVerticesOld = obj.freeVertices;
+ obj.freeVertices = getFreeDofs(obj.fes);
+
+ nonInvertedSmootherMatrix = A;
+ obj.inverseDiagonal{L} = full(diag(A)).^(-1);
+
+ if L >= 2
+ obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
+ obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
+ obj.inverseDiagonal{L} = obj.inverseDiagonal{L}(obj.changedPatches{L});
+ end
+ end
+
+ function Cx = smooth(obj, x, k)
+ idx = obj.changedPatches{k};
+ Cx = zeros(size(x));
+ Cx(idx,:) = obj.inverseDiagonal{k} .* x(idx,:);
+ end
+
+ function Px = prolongate(obj, x, k)
+ Px = obj.intergridMatrix{k} * x;
+ end
+
+ function Px = restrict(obj, x, k)
+ Px = obj.intergridMatrix{k}' * x;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/spaces/@FeSpace/FeSpace.m b/lib/spaces/@FeSpace/FeSpace.m
index 357fb167d8f0ab37b57a72b631843a5fcd93e7ff..39413fd4187254266dcd80e86aa98cad588a2692 100644
--- a/lib/spaces/@FeSpace/FeSpace.m
+++ b/lib/spaces/@FeSpace/FeSpace.m
@@ -84,6 +84,10 @@ classdef FeSpace < handle
end
freeDofs = obj.freeDofs;
end
+
+ dofs = createVertexDofs(obj, idx)
+ dofs = createInnerEdgeDofs(obj, idx)
+ dofs = createInnerElementDofs(obj, idx)
end
methods (Access=protected)
diff --git a/lib/spaces/@FeSpace/assembleDofs.m b/lib/spaces/@FeSpace/assembleDofs.m
index 37a867db465801e8a0fc2e76234800122960a461..ca5ad9ecc4c5d39f0fe6e3e41dbc63340888c717 100644
--- a/lib/spaces/@FeSpace/assembleDofs.m
+++ b/lib/spaces/@FeSpace/assembleDofs.m
@@ -10,41 +10,34 @@ nDofs = dofCoupling * [3;3;1];
nEdgeDofs = dofCoupling * [2;1;0];
dofs.element2Dofs = zeros(nDofs, mesh.nElements);
dofs.edge2Dofs = zeros(nEdgeDofs, mesh.nEdges);
-n = 0; % current number of local ...
-k = 0; % ... local edge ...
-N = 0; % ... and global Dofs
+
+dofs.nDofs = dofCoupling * [mesh.nCoordinates; mesh.nEdges; mesh.nElements];
%% vertex Dofs
if dofCoupling(1) ~= 0
- dofs.element2Dofs(1:3,:) = mesh.elements;
- dofs.edge2Dofs(1:2,:) = mesh.edges;
- k = k + 2;
- n = n + 3;
- N = N + mesh.nCoordinates;
+ dofNumbers = createVertexDofs(obj);
+ dofs.element2Dofs(1:3,:) = dofNumbers(mesh.elements);
+ dofs.edge2Dofs(1:2,:) = dofNumbers(mesh.edges);
end
%% edge dofs
if dofCoupling(2) ~= 0
nInnerEdgeDofs = dofCoupling(2);
- innerEdgeDofs = N + reshape(1:(nInnerEdgeDofs*mesh.nEdges), nInnerEdgeDofs, mesh.nEdges);
- for i = 1:3
- idx = (n+1):(n+nInnerEdgeDofs);
- dofs.element2Dofs(idx,:) = innerEdgeDofs(:,mesh.element2edges(i,:));
- dofs.element2Dofs(idx,mesh.flipEdges(i,:)) = flipud(dofs.element2Dofs(idx,mesh.flipEdges(i,:)));
- n = n + nInnerEdgeDofs;
+ dofNumbers = createInnerEdgeDofs(obj);
+ idx = dofCoupling(1) * 3;
+ for k = 1:3
+ idx = idx(end) + (1:nInnerEdgeDofs);
+ dofs.element2Dofs(idx,:) = dofNumbers(:,mesh.element2edges(k,:));
+ dofs.element2Dofs(idx,mesh.flipEdges(k,:)) = flipud(dofs.element2Dofs(idx,mesh.flipEdges(k,:)));
end
- dofs.edge2Dofs((k+1):(k+nInnerEdgeDofs),:) = innerEdgeDofs;
- N = N + nInnerEdgeDofs*mesh.nEdges;
+ idx = dofCoupling(1) * 2 + (1:nInnerEdgeDofs);
+ dofs.edge2Dofs(idx,:) = dofNumbers;
end
%% inner dofs
if dofCoupling(3) ~= 0
- nInnerDofsPerElement = dofCoupling(3);
- dofs.element2Dofs((n+1):(n+nInnerDofsPerElement),:) ...
- = N + reshape(1:(nInnerDofsPerElement*mesh.nElements), nInnerDofsPerElement, mesh.nElements);
- N = N + nInnerDofsPerElement*mesh.nElements;
+ idx = dofCoupling(1:2) * [3;3] + (1:dofCoupling(3));
+ dofs.element2Dofs(idx,:) = createInnerElementDofs(obj);
end
-dofs.nDofs = N;
-
end
diff --git a/lib/spaces/@FeSpace/createInnerEdgeDofs.m b/lib/spaces/@FeSpace/createInnerEdgeDofs.m
new file mode 100644
index 0000000000000000000000000000000000000000..35a8ddca1f09d9d35153ab6cad1bc552b055bedf
--- /dev/null
+++ b/lib/spaces/@FeSpace/createInnerEdgeDofs.m
@@ -0,0 +1,38 @@
+% createInnerEdgeDofs Generate global DOF-numbers associated with inner edge
+% dofs.
+%
+% dofs = createInnerEdgeDofs(obj) generates dofs for all edges.
+%
+% dofs = createInnerEdgeDofs(obj, idx) generates dofs only for edges given by
+% index vector idx.
+
+function dofs = createInnerEdgeDofs(obj, idx)
+
+arguments
+ obj
+ idx (1,:) {mustBeIndexVector} = ':'
+end
+
+N = obj.mesh.nEdges;
+nLocalDofs = getDofConnectivity(obj.finiteElement);
+
+if nLocalDofs(2) == 0
+ dofs = [];
+ return
+end
+
+% dof numbers are simply given by vertex numbers
+if ischar(idx) && isequal(idx, ':')
+ dofs = 1:N;
+elseif islogical(idx)
+ assert(isequal(size(idx), [1, N]))
+ dofs = find(idx);
+else
+ assert(all(0 < idx) && all(idx <= N))
+ dofs = idx;
+end
+
+offset = obj.mesh.nCoordinates * nLocalDofs(1);
+dofs = offset + (nLocalDofs(2)*(dofs-1) + (1:nLocalDofs(2))');
+
+end
diff --git a/lib/spaces/@FeSpace/createInnerElementDofs.m b/lib/spaces/@FeSpace/createInnerElementDofs.m
new file mode 100644
index 0000000000000000000000000000000000000000..93a41d41af9838079c37dde598ba327933a2a0be
--- /dev/null
+++ b/lib/spaces/@FeSpace/createInnerElementDofs.m
@@ -0,0 +1,38 @@
+% createInnerElementDofs Generate global DOF-numbers associated with inner
+% element dofs.
+%
+% dofs = createInnerElementDofs(obj) generates dofs for all elements.
+%
+% dofs = createInnerElementDofs(obj, idx) generates dofs only for elements
+% given by index vector idx.
+
+function dofs = createInnerElementDofs(obj, idx)
+
+arguments
+ obj
+ idx (1,:) {mustBeIndexVector} = ':'
+end
+
+N = obj.mesh.nElements;
+nLocalDofs = getDofConnectivity(obj.finiteElement);
+
+if nLocalDofs(3) == 0
+ dofs = [];
+ return
+end
+
+% dof numbers are simply given by vertex numbers
+if ischar(idx) && isequal(idx, ':')
+ dofs = 1:N;
+elseif islogical(idx)
+ assert(isequal(size(idx), [1, N]))
+ dofs = find(idx);
+else
+ assert(all(0 < idx) && all(idx <= N))
+ dofs = idx;
+end
+
+offset = nLocalDofs(1:2) * [obj.mesh.nCoordinates; obj.mesh.nEdges];
+dofs = offset + (nLocalDofs(3)*(dofs-1) + (1:nLocalDofs(3))');
+
+end
diff --git a/lib/spaces/@FeSpace/createVertexDofs.m b/lib/spaces/@FeSpace/createVertexDofs.m
new file mode 100644
index 0000000000000000000000000000000000000000..d8985325cdf623718cbcc192cf546a19f202168a
--- /dev/null
+++ b/lib/spaces/@FeSpace/createVertexDofs.m
@@ -0,0 +1,34 @@
+% createVertexDofs Generate global DOF-numbers associated with vertex dofs.
+%
+% dofs = createVertexDofs(obj) generates dofs for all vertices.
+%
+% dofs = createVertexDofs(obj, idx) generates dofs only for vertices given by
+% index vector idx.
+
+function dofs = createVertexDofs(obj, idx)
+
+arguments
+ obj
+ idx (1,:) {mustBeIndexVector} = ':'
+end
+
+N = obj.mesh.nCoordinates;
+nLocalDofs = getDofConnectivity(obj.finiteElement);
+
+if nLocalDofs(1) == 0
+ dofs = [];
+ return
+end
+
+% dof numbers are simply given by vertex numbers
+if ischar(idx) && isequal(idx, ':')
+ dofs = 1:N;
+elseif islogical(idx)
+ assert(isequal(size(idx), [1, N]))
+ dofs = find(idx);
+else
+ assert(all(0 < idx) && all(idx <= N))
+ dofs = idx;
+end
+
+end
diff --git a/lib/spaces/interpolation/LinearInterpolation.m b/lib/spaces/interpolation/LinearInterpolation.m
index d5d1e7db815f156d168e6d17cb4d3fe2a0844ab2..02409a6f344140117475d2c6ff1365985232da3f 100644
--- a/lib/spaces/interpolation/LinearInterpolation.m
+++ b/lib/spaces/interpolation/LinearInterpolation.m
@@ -1,7 +1,5 @@
% LinearInterpolation (subclass of Evaluable) Provide piecewise linear
-% interpolant of function. This is used mainly internally for prolongation.
-%
-% See also: LinearFeProlongation
+% interpolant of function. This can be used for approximate prolongation.
classdef LinearInterpolation < Evaluable
%% properties
diff --git a/lib/spaces/interpolation/LoFeProlongation.m b/lib/spaces/prolongation/LoMeshProlongation.m
similarity index 84%
rename from lib/spaces/interpolation/LoFeProlongation.m
rename to lib/spaces/prolongation/LoMeshProlongation.m
index 2c2603e8cf42a7371dabb010b93874d389047f30..b6ebbda297ba295c74b35676da1b6077d78484b9 100644
--- a/lib/spaces/interpolation/LoFeProlongation.m
+++ b/lib/spaces/prolongation/LoMeshProlongation.m
@@ -1,29 +1,30 @@
-% LoFeProlongation (subclass of Prolongation) Prolongate lowest order L^2/H^1
+% LoMeshProlongation (subclass of MeshProlongation) Prolongate lowest order L^2/H^1
% conforming finite element function to refined mesh.
%
-% P = LoFeProlongation(fes) returns a handle to the prolongation object
+% P = LoMeshProlongation(fes) returns a handle to the prolongation object
% associated to the finite element space fes. The prolongation matrix
% P.matrix is set automatically at mesh refinement.
%
-% prolongate(P, u) returns the prolongated data of FeFunction u.
+% See also Prolongation
-classdef LoFeProlongation < Prolongation
+classdef LoMeshProlongation < MeshProlongation
%% properties
properties (Access=protected)
feType
end
+
%% methods
methods (Access=public)
- function obj = LoFeProlongation(fes)
- obj = obj@Prolongation(fes);
+ function obj = LoMeshProlongation(fes)
+ obj = obj@MeshProlongation(fes);
switch class(fes.finiteElement)
case 'LowestOrderL2Fe'
obj.feType = 'L2';
case 'LowestOrderH1Fe'
obj.feType = 'H1';
otherwise
- eid = 'LoFeProlongation:wrongFeType';
- msg = 'LoFeProlongation needs a lowest order L2 or H1 finite element space.';
+ eid = 'LoMeshProlongation:wrongFeType';
+ msg = 'LoMeshProlongation needs a lowest order L2 or H1 finite element space.';
throwAsCaller(MException(eid, msg));
end
end
@@ -87,5 +88,9 @@ classdef LoFeProlongation < Prolongation
sum(data.nInnerNodes.*data.nRefinedElements);
obj.matrix = sparse(I, J, V, nNewDofs, dofs.nDofs);
end
+
+ % no-op override to avoid wrong call to superclass method
+ function connectDofs(~, ~, ~)
+ end
end
end
diff --git a/lib/spaces/interpolation/FeProlongation.m b/lib/spaces/prolongation/MeshProlongation.m
similarity index 84%
rename from lib/spaces/interpolation/FeProlongation.m
rename to lib/spaces/prolongation/MeshProlongation.m
index 4dc5d6b97ac4e223403a4794c378ad118f6e3a8a..efc1c76e4366696b9540258fe17c1d06a48cd296 100644
--- a/lib/spaces/interpolation/FeProlongation.m
+++ b/lib/spaces/prolongation/MeshProlongation.m
@@ -1,22 +1,31 @@
-% FeProlongation (subclass of Prolongation) Provides prolongation operator for
+% MeshProlongation (subclass of Prolongation) Provides prolongation operator for
% general FeFunction to a refined mesh.
%
-% P = FeProlongation(fes) returns a handle to the prolongation object
+% P = MeshProlongation(fes) returns a handle to the prolongation object
% associated to the finite element space fes. The prolongation matrix
% P.matrix is set automatically at mesh refinement.
%
-% prolongate(P, u) returns the prolongated data of FeFunction u.
+% See also: Prolongation
-classdef FeProlongation < Prolongation
+classdef MeshProlongation < Prolongation
%% properties
properties (Access=private)
postRefineListener
end
+
+ properties (Access=protected)
+ fes FeSpace
+ listenerHandle
+ end
%% methods
methods (Access=public)
- function obj = FeProlongation(fes)
- obj = obj@Prolongation(fes);
+ function obj = MeshProlongation(fes)
+ obj = obj@Prolongation();
+ obj.fes = fes;
+
+ mesh = fes.mesh;
+ obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.setupMatrix);
obj.postRefineListener = fes.mesh.listener('RefineCompleted', @obj.connectDofs);
assert(isa(obj.fes.finiteElement, 'NodalFiniteElement'), ...
@@ -26,10 +35,8 @@ classdef FeProlongation < Prolongation
methods (Access=protected)
function setupMatrix(obj, ~, data)
- % general idea: compute *all* dofs for each new element and store
- % them consecutively
- % those are connected to the actual new dofs when their numbering
- % per element is available (-> obj.connectDofs)
+ % general idea: compute *all* dofs for each new element and store them consecutively
+ % those are connected to the actual new dofs when their numbering per element is available (-> obj.connectDofs)
% de-allocate for memory-efficiency
obj.matrix = [];
@@ -61,8 +68,7 @@ classdef FeProlongation < Prolongation
for k = find(data.nRefinedElements ~= 0)'
unitTriangle = getBisectedUnitTriangle(class(data.bisection{k}));
newDofLocations = getNewLocationsElementwise(unitTriangle, dofLocations);
- newDofWeights = divideElementwise(...
- reshape(evalShapeFunctions(fe, newDofLocations), nLocalDofs, nLocalDofs, []), localMatrix);
+ newDofWeights = divideElementwise(evalShapeFunctions(fe, newDofLocations), localMatrix);
nNewLocs = newDofLocations.nNodes;
elems = getRefinedElementIdx(data, k);
@@ -91,7 +97,7 @@ end
%% local functions
function mesh = getBisectedUnitTriangle(bisecMethod)
- mesh = Mesh.loadFromGeometry('unittriangle');
+ mesh = Mesh.unitTriangle();
% bisect unit triangle according to given bisection rule
% TODO: this should be handled by bisection method itself?
@@ -123,10 +129,8 @@ function elementwiseLocations = getNewLocationsElementwise(mesh, locations)
end
function C = divideElementwise(A, B)
- % TODO: from 2022a on, this can be replaced by pagemrdivide(A, B)
n = size(A,1);
- C = B' \ reshape(permute(A, [2,1,3]), n, []);
- C = reshape(permute(reshape(C, n, n, []), [2,1,3]), n, []);
+ C = reshape(pagemrdivide(reshape(A, n, n, []), B), n, []);
end
function idx = getConsecutiveIndices(parents, nIndices)
diff --git a/lib/spaces/interpolation/Prolongation.m b/lib/spaces/prolongation/Prolongation.m
similarity index 55%
rename from lib/spaces/interpolation/Prolongation.m
rename to lib/spaces/prolongation/Prolongation.m
index 97adb20868910b445d6a1884181623ff17abeab0..f6ca211341438e536a09b9776e0c4121647c2912 100644
--- a/lib/spaces/interpolation/Prolongation.m
+++ b/lib/spaces/prolongation/Prolongation.m
@@ -1,5 +1,9 @@
% Prolongation (abstract handle class) Interface for prolongation of FeFunctions
% from coarse to fine mesh.
+%
+% prolongate(P, u) returns the prolongated data of FeFunction u.
+%
+% restrict(P, u) returns the restricted data of FeFunction u.
classdef Prolongation < handle
%% properties
@@ -7,18 +11,11 @@ classdef Prolongation < handle
matrix {mustBeSparse}
end
- properties (Access=protected)
- fes FeSpace
- listenerHandle
- end
%% methods
methods (Access=public)
function obj = Prolongation(fes)
- obj.fes = fes;
- mesh = fes.mesh;
obj.matrix = [];
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.setupMatrix);
end
function newFeData = prolongate(obj, u)
@@ -32,10 +29,29 @@ classdef Prolongation < handle
end
newFeData = obj.matrix * u.data';
end
+
+ function newFeData = restrict(obj, u)
+ arguments
+ obj
+ u FeFunction
+ end
+
+ if isempty(obj.matrix)
+ error('Prolongation matrix not set. Are you sure there was mesh refinement?')
+ end
+ newFeData = (u.data * obj.matrix)';
+ end
end
- methods (Abstract, Access=protected)
- setupMatrix(obj, src, event)
+ methods (Static, Access=public)
+ function P = chooseFor(fes)
+ if isa(fes.finiteElement, 'LowestOrderH1Fe') ...
+ || isa(fes.finiteElement, 'LowestOrderL2Fe')
+ P = LoMeshProlongation(fes);
+ else
+ P = MeshProlongation(fes);
+ end
+ end
end
end
diff --git a/lib/spaces/prolongation/SpaceProlongation.m b/lib/spaces/prolongation/SpaceProlongation.m
new file mode 100644
index 0000000000000000000000000000000000000000..83d690006d817d2a48c865941be819e44386eaa4
--- /dev/null
+++ b/lib/spaces/prolongation/SpaceProlongation.m
@@ -0,0 +1,57 @@
+% SpaceProlongation (subclass of Prolongation) Provides prolongation operator
+% from S^p to S^q with q > p
+%
+% P = SpaceProlongation(sourceFes, targetFes) returns a handle to the
+% prolongation object for the prolongation sourceFes -> targetFes. The
+% prolongation matrix P.matrix is set up at construction.
+%
+% See also: Prolongation
+
+classdef SpaceProlongation < Prolongation
+ %% methods
+ methods (Access=public)
+ function obj = SpaceProlongation(sourceFes, targetFes)
+ unitTriangleInclusionMatrix = SpaceProlongation.getUnitTriangleInclusionMatrix(sourceFes.finiteElement, targetFes.finiteElement);
+
+ % Copy the matrix from the unit triangle for each element
+ mat = repmat(unitTriangleInclusionMatrix(:), [1, sourceFes.mesh.nElements]);
+
+ % Create index sets for the accumarray
+ fromDofs = getDofs(sourceFes);
+ toDofs = getDofs(targetFes);
+ I = repmat(fromDofs.element2Dofs, size(toDofs.element2Dofs, 1), 1);
+ J = repelem(toDofs.element2Dofs, size(fromDofs.element2Dofs, 1), 1);
+
+ % assemble matrix and account for multiply computed dofs
+ obj.matrix = sparse(I(:), J(:), mat(:), fromDofs.nDofs, toDofs.nDofs);
+ dofCounts = accumarray(toDofs.element2Dofs(:), 1, [toDofs.nDofs, 1])';
+ obj.matrix = obj.matrix ./ dofCounts;
+ end
+ end
+
+ methods (Static, Access=protected)
+ function matrix = getUnitTriangleInclusionMatrix(sourceFE, targetFE)
+ unittriangle = Mesh.unitTriangle();
+ fromFes = FeSpace(unittriangle, sourceFE);
+ toFes = FeSpace(unittriangle, targetFE);
+ Vertices = 1:getDofs(fromFes).nDofs;
+ matrix = SpaceProlongation.interpolateData(eye(length(Vertices)), fromFes, toFes)';
+
+ % Get local to global map in unit triangle
+ matrix = matrix(getDofs(fromFes).element2Dofs, getDofs(toFes).element2Dofs);
+ end
+
+ function interpolatedData = interpolateData(data, fromFes, toFes)
+ Dofs = 1:getDofs(toFes).nDofs;
+ nComponents = size(data, 2);
+ interpolatedData = zeros(numel(Dofs), nComponents);
+ feFunctionWrapper = FeFunction(fromFes);
+ for k = 1:nComponents
+ feFunctionWrapper.setData(data(:,k));
+ wholeData = nodalInterpolation(feFunctionWrapper, toFes);
+ interpolatedData(:,k) = wholeData(Dofs);
+ end
+ end
+ end
+end
+
diff --git a/lib/storage/@LevelData/LevelData.m b/lib/storage/@LevelData/LevelData.m
new file mode 100644
index 0000000000000000000000000000000000000000..5b4955d15d2da5300edffc333fa6568f8d35ca78
--- /dev/null
+++ b/lib/storage/@LevelData/LevelData.m
@@ -0,0 +1,299 @@
+classdef LevelData < handle & matlab.mixin.CustomDisplay
+%%LEVELDATA class representing results from level-oriented computations
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ properties
+ % map for general metadata
+ metaData (1,1) dictionary
+ % Root folder for file storage
+ root
+ % Structure array storing string representations of variables
+ dictionary
+ % Array of level data
+ level (1,:) struct
+ % Dictionary of data types
+ type (1,1) dictionary
+ end
+
+ properties (Access=private)
+ % Dictionary of categories for cataloguing and plotting
+ category (1,1) dictionary
+ end
+
+ properties (Dependent)
+ % Path to folder for file storage
+ foldername (1,:) char
+ % Name of file for storage
+ filename (1,:) char
+ % Cell array of variable names
+ label (1,:) cell
+ % Number of levels
+ nLevel (1,1) int32
+ % Number of variables
+ nVariable (1,1) int32
+ % Number of scalar variables
+ nScalarVariable (1,1) int32
+ % Cell array of names of scalar variables
+ scalarVariable (1,:) cell
+ % Number of absolute value variables
+ nAbsoluteVariable (1,1) int32
+ % Number of time variables
+ nTimeVariable (1,1) int32
+ % Boolean indicating whether data has been stored for a single level
+ isInitialRun (1,1) logical
+ % Cell array of names of absolute value variables
+ absoluteVariable (1,:) cell
+ % Cell array of names of timing variables
+ timeVariable (1,:) cell
+ end
+
+ properties (Dependent, Access=private)
+ % Boolean vector determining whether variable is a scalar value
+ isScalar (1,:) logical
+ end
+
+ methods
+ %% CONSTRUCTOR
+ function obj = LevelData(rootpath)
+ %%LEVELDATA creates an object for storing level-oriented data,
+ %file storage is initialised in the specified optional path
+ % leveldata = LEVELDATA(rootpath)
+
+ arguments
+ rootpath {mustBeTextScalar} = 'results'
+ end
+ obj.root = rootpath;
+
+ % TODO: set output of hostname to string and remove cast
+ obj.metaData = dictionary(...
+ "problem", "problem", ...
+ "domain", "domain", ...
+ "method", "method", ...
+ "identifier", "main", ...
+ "hostname", string(getHostname()), ...
+ "timestamp", string(datetime('now', 'Format', 'yyyy-MM-dd_HH:mm:ss')));
+
+ obj.category = dictionary();
+
+ % Initialise dictionary with some default values
+ obj.dictionary = getDefaultDictionary();
+ end
+
+ %% SET GLOBAL VARIABLES
+ function set.root(obj, path)
+ arguments
+ obj
+ path {mustBeTextScalar}
+ end
+ obj.root = path;
+ end
+
+ %% GET GLOBAL VARIABLES
+ function label = get.label(obj)
+ % List of names of all variables
+ if obj.nLevel > 0
+ label = fieldnames(obj.level(1));
+ else
+ label = {};
+ end
+ end
+
+ function nLevel = get.nLevel(obj)
+ nLevel = length(obj.level);
+ end
+
+ function nVariable = get.nVariable(obj)
+ nVariable = length(obj.label);
+ end
+
+ function path = get.foldername(obj)
+ % Name of folder for file storage
+ rootpath = obj.root;
+ if numel(rootpath) > 0 && ~strcmp(rootpath(end), '/')
+ rootpath = [rootpath, '/'];
+ end
+ path = rootpath + strjoin([obj.metaData("problem"), obj.metaData("domain"), obj.metaData("method")], '_');
+ end
+
+ function file = get.filename(obj)
+ % Name of file for file storage
+ file = [obj.metaData("identifier")];
+ end
+
+ function bool = get.isScalar(obj)
+ % Creates boolean vector determining whether a variable
+ % contains a scalar value per level only
+ bool = [obj.type(obj.label).isScalar]';
+ end
+
+ function variables = get.scalarVariable(obj)
+ % Returns list of names of all scalar variables
+ variables = obj.label(obj.isScalar);
+ end
+
+ function nVariable = get.nScalarVariable(obj)
+ nVariable = nnz(obj.isScalar);
+ end
+
+ function variables = get.timeVariable(obj)
+ % Returns list of names of all time variables
+ idx = (obj.category(obj.label) == DataCategory.TIME);
+ variables = obj.label(idx);
+ end
+
+ function nVariable = get.nTimeVariable(obj)
+ nVariable = nnz(obj.category.values() == DataCategory.TIME);
+ end
+
+ function variables = get.absoluteVariable(obj)
+ % Returns list of names of all absolute variables
+ idx = (obj.category(obj.label) == DataCategory.ABSOLUTE);
+ variables = obj.label(idx);
+ end
+
+ function nVariable = get.nAbsoluteVariable(obj)
+ nVariable = nnz(obj.category.values() == DataCategory.ABSOLUTE);
+ end
+
+ function bool = get.isInitialRun(obj)
+ % True if only one level is stored
+ bool = (obj.nLevel <= 1);
+ end
+
+
+ %% READ LEVEL DATA
+ data = get(obj, jLevel, variableName)
+
+ %% MODIFY LEVEL DATA
+ set(obj, jLevel, variableName, value)
+ setTime(obj, jLevel, variableName, value)
+ setAbsolute(obj, jLevel, variableName, value)
+
+ function append(obj, varargin)
+ %%APPEND simplified addition of specified data for the
+ %subsequent level, the arguments must be specified as
+ %pairs of variable names and data (numeric or string /
+ %characters)
+ % APPEND(obj, variableName, value, ...)
+
+ arguments
+ obj
+ end
+
+ arguments (Repeating)
+ varargin
+ end
+
+ assert(mod(length(varargin), 2) == 0, ...
+ 'Invalid number of arguments');
+ nlvl = obj.nLevel + 1;
+ obj.set(nlvl, varargin{:});
+ end
+
+ function remove(obj, variableName)
+ %%REMOVE removes the specified list of variables
+ % REMOVE(obj, variableName, ...)
+
+ arguments
+ obj
+ end
+
+ arguments (Repeating)
+ variableName {mustBeTextScalar}
+ end
+
+ obj.level = rmfield(obj.level, variableName);
+ obj.type(variableName) = [];
+ end
+
+ function removeLevel(obj, jLevel)
+ %%REMOVELEVEL removes all data for the specified list jLevel
+ %of level numbers
+ % REMOVELEVEL(obj, jLevel)
+ for k = 1:length(jLevel)
+ obj.level(jLevel(k)) = [];
+ end
+ end
+
+ function removeNonscalar(obj)
+ %%REMOVENONSCALAR removes all non-scalar data completely;
+ %e.g., to reduce storage size
+ % REMOVENONSCALAR(obj)
+ variableNonScalar = setdiff(obj.label, obj.scalarVariable);
+ obj.remove(variableNonScalar{:});
+ end
+
+ %% OUTPUT LEVEL DATA
+ printHeader(obj)
+ printLevel(obj, jLevel)
+ ax = plot(obj, xVariable, yVariable)
+ ax = plotTime(obj, xVariable, yVariable)
+ ax = plotAbsolute(obj, xVariable, yVariable)
+ ax = plotTriangulation(obj, jLevel)
+
+ %% EXPORT LEVEL DATA
+ saveToFile(obj, folder, file)
+ saveToTable(obj, separator)
+ plotToFile(obj, xVariable, yVariable)
+ plotTriangulationToFile(obj, jLevel)
+ end
+
+ methods (Access=protected)
+ %% CUSTOM DISPLAY
+ propgrp = getPropertyGroups(obj)
+ end
+
+ %% AUXILIARY FUNCTIONS
+ methods (Access = private)
+ function idx = getIndex(obj, variableName)
+ [~, idx] = ismember(variableName, obj.label);
+ end
+
+ function width = getWidth(~, type)
+ minimalWidth = 8;
+ width = num2str(max(type.printWidth, minimalWidth));
+ end
+
+ function spec = getHeaderSpecifier(obj, separator)
+ % Creates formatting string for the header of the output to command line
+ spec = cell(1, obj.nScalarVariable);
+ for j = 1:obj.nScalarVariable
+ t = obj.type(obj.scalarVariable{j});
+ spec{j} = assembleSpecifier(obj.getWidth(t), 's');
+ end
+ spec = strjoin(spec, separator) + "\n";
+ end
+
+ function spec = getFormatSpecifier(obj, separator)
+ % Creates formatting string for printing to command line
+ spec = cell(1, obj.nScalarVariable);
+ for j = 1:obj.nScalarVariable
+ t = obj.type(obj.scalarVariable{j});
+ spec{j} = assembleSpecifier(obj.getWidth(t), t.formatSpec);
+ end
+ spec = strjoin(spec, separator) + "\n";
+ end
+
+ ax = plotLevel(obj, plotFunction, xVariable, variableName)
+ plotLevelTriangulation(obj, ax, jLevel)
+ end
+end
+
+function spec = assembleSpecifier(width, format)
+ spec = ['%', num2str(width), format];
+end
diff --git a/lib/storage/@LevelData/get.m b/lib/storage/@LevelData/get.m
new file mode 100644
index 0000000000000000000000000000000000000000..c130089623ad775a2c9d93b97a96c12135fc96ff
--- /dev/null
+++ b/lib/storage/@LevelData/get.m
@@ -0,0 +1,74 @@
+function data = get(obj, jLevel, variableName)
+%%GET extracts data from this LevelData object on a given list jLevel of
+%level numbers. One or more variables to be returned can be specified by their
+%names.
+% data = GET(obj, jLevel, variableName, ...)
+% data = GET(obj, ':', variableName, ...)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jLevel {mustBeIndexVector} = ':'
+ end
+
+ arguments (Repeating)
+ variableName {mustBeTextScalar}
+ end
+
+ if strcmp(jLevel, ':')
+ jLevel = 1:obj.nLevel;
+ end
+
+ % Initialise return variable
+ data = nan(length(jLevel), length(variableName));
+
+ % filter non-existing variables
+ existingVariables = ismember(variableName, obj.label);
+ if any(~existingVariables)
+ warning(['Variable(s) ', strjoin(variableName(~existingVariables), ', '), ' not found']);
+ end
+ variableName = variableName(existingVariables);
+
+ % Iterate over all variables
+ for jVariable = 1:length(variableName)
+ name = variableName{jVariable};
+
+ % Determine index of current variable
+ idx = obj.getIndex(name);
+ if obj.isScalar(idx)
+ % Extract scalar variables for each level
+ value = [obj.level.(name)];
+ else
+ % Extract non-scalar variables only in case of a single level
+ if length(jLevel) > 1
+ warning('Export of level-oriented scalar variables only');
+ value = [];
+ else
+ data(1,jVariable) = obj.level(jLevel).(name);
+ return
+ end
+ end
+
+ % Save extracted data to return variable
+ data(:,jVariable) = value;
+ % Post-process character arrays
+ if obj.type(name).rawType == RawType.TEXT
+ data = mat2cell(char(data), ones(length(data), 1));
+ end
+ end
+end
diff --git a/lib/storage/@LevelData/getPropertyGroups.m b/lib/storage/@LevelData/getPropertyGroups.m
new file mode 100644
index 0000000000000000000000000000000000000000..3063a08ad7e4e79ddc2c4a74372b35879995a836
--- /dev/null
+++ b/lib/storage/@LevelData/getPropertyGroups.m
@@ -0,0 +1,24 @@
+function propgrp = getPropertyGroups(obj)
+ if ~isscalar(obj)
+ propgrp = getPropertyGroups@matlab.mixin.CustomDisplay(obj);
+ else
+ % some rudimentary statistics about variable counts
+ statList = ["nLevel", "nVariable", "nScalarVariable", "nAbsoluteVariable", "nTimeVariable"];
+ statTitle = "Variable Counts";
+ statGrp = matlab.mixin.util.PropertyGroup(statList, statTitle);
+
+ % metadata as a struct
+ metaDataList = struct();
+ for k = obj.metaData.keys'
+ metaDataList.(k{1}) = obj.metaData(k{1});
+ end
+ metaDataTitle = "Metadata";
+ metaDataGrp = matlab.mixin.util.PropertyGroup(metaDataList, metaDataTitle);
+
+ % storage related properties
+ storeList = ["root", "foldername", "filename"];
+ storeTitle = "Storage details";
+ storeGrp = matlab.mixin.util.PropertyGroup(storeList, storeTitle);
+ propgrp = [statGrp, metaDataGrp, storeGrp];
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/plot.m b/lib/storage/@LevelData/plot.m
new file mode 100644
index 0000000000000000000000000000000000000000..1151963f6c7313b49169568fbd5fa489971ca05b
--- /dev/null
+++ b/lib/storage/@LevelData/plot.m
@@ -0,0 +1,38 @@
+function ax = plot(obj, xVariable, yVariable)
+%%PLOT plots various scalar variables (specified in yVariable) with respect
+%to the variable xVariable, returns handle to axis object
+% ax = PLOT(obj, xVariable, yVariable, ...)
+%
+% See also LevelData/plotTime, LevelData/plotAbsolute
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ xVariable {mustBeTextScalar}
+ end
+
+ arguments (Repeating)
+ yVariable {mustBeTextScalar}
+ end
+
+ if isempty(yVariable)
+ yVariable = setdiff(obj.label, xVariable);
+ end
+
+ ax = obj.plotLevel(DataCategory.ERROR, xVariable, yVariable);
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/plotAbsolute.m b/lib/storage/@LevelData/plotAbsolute.m
new file mode 100644
index 0000000000000000000000000000000000000000..931353fe3b137c3be13c32a0fe90d7ceeac0eaa0
--- /dev/null
+++ b/lib/storage/@LevelData/plotAbsolute.m
@@ -0,0 +1,39 @@
+function ax = plotAbsolute(obj, xVariable, yVariable)
+%%PLOTABSOLUTE plots various scalar absolute variables (specified in
+%yVariable) with respect to the variable xVariable, returns handle to axis
+%object
+% ax = PLOTABSOLUTE(obj, xVariable, yVariable, ...)
+%
+% See also LevelData/plot, LevelData/plotTime
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ xVariable {mustBeTextScalar}
+ end
+
+ arguments (Repeating)
+ yVariable {mustBeTextScalar}
+ end
+
+ if isempty(yVariable)
+ yVariable = setdiff(obj.label, xVariable);
+ end
+
+ ax = obj.plotLevel(DataCategory.ABSOLUTE, xVariable, yVariable);
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/plotLevel.m b/lib/storage/@LevelData/plotLevel.m
new file mode 100644
index 0000000000000000000000000000000000000000..455d8290b990e4e7597c84b78ffebf60232c673a
--- /dev/null
+++ b/lib/storage/@LevelData/plotLevel.m
@@ -0,0 +1,74 @@
+function ax = plotLevel(obj, category, xVariable, yVariable)
+%%PLOTLEVEL auxiliary private function for creation of plots for usage in
+%LevelData/plot, LevelData/plotTime, and LevelData/plotAbsolute, uses the
+%specified DataCategory for the plot of variables in yVariable with
+%respect to xVariable, returns handle to axis object
+% ax = PLOTLEVEL(obj, category, xVariable, yVariable, ...)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ category DataCategory
+ xVariable {mustBeTextScalar}
+ yVariable (1,:) cell
+ end
+
+ % Extract variables with correct category, type, and shape for y-axis
+ idx = (obj.category(yVariable) == category) ...
+ & ([obj.type(yVariable).rawType] == RawType.FLOAT) ...
+ & [obj.type(yVariable).isScalar];
+ yVariable = yVariable(idx);
+
+ % Create handle to currently active axis object
+ ax = gca;
+
+ % Extract variable for x-axis
+ xValue = obj.get(1:obj.nLevel, xVariable);
+
+ % Create unified list for formatting plot lines
+ [COLOURORDER, MARKER] = getPlotStyle();
+
+ % Iterate over given variables
+ for j = 1:length(yVariable)
+ % Extract value for y-axis
+ yValue = obj.get(1:obj.nLevel, yVariable{j});
+
+ % Extract label for legend from dictionary
+ if obj.dictionary.isKey(yVariable{j})
+ variableLabel = obj.dictionary(yVariable{j});
+ else
+ variableLabel = yVariable{j};
+ end
+
+ % Create plot
+ category.plotFunction( ...
+ ax, xValue, yValue, '-', ...
+ 'Marker', MARKER{mod(j, length(MARKER))}, ...
+ 'Color', COLOURORDER(j, :), ...
+ 'LineWidth', 1.5, ...
+ 'DisplayName', variableLabel);
+ % Add new line into the current figure when calling plotConvergence again
+ hold(ax, 'on');
+ end
+
+ % Set plot appearance
+ title(ax, category.title);
+ xlabel(ax, xVariable);
+ ylabel(ax, category.yLabel);
+ configureLegend(ax, category.legendLocation);
+end
diff --git a/lib/storage/@LevelData/plotTime.m b/lib/storage/@LevelData/plotTime.m
new file mode 100644
index 0000000000000000000000000000000000000000..3b7693f5f684259fe1465bff7b0ccdff03de12fc
--- /dev/null
+++ b/lib/storage/@LevelData/plotTime.m
@@ -0,0 +1,38 @@
+function ax = plotTime(obj, xVariable, yVariable)
+%%PLOTTIME plots various scalar time variables (specified in yVariable) with
+%respect to the variable xVariable, returns handle to axis object
+% ax = PLOTTIME(obj, xVariable, yVariable, ...)
+%
+% See also LevelData/plot, LevelData/plotAbsolute
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ xVariable {mustBeTextScalar}
+ end
+
+ arguments (Repeating)
+ yVariable {mustBeTextScalar}
+ end
+
+ if isempty(yVariable)
+ yVariable = setdiff(obj.label, xVariable);
+ end
+
+ ax = obj.plotLevel(DataCategory.TIME, xVariable, yVariable);
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/plotToFile.m b/lib/storage/@LevelData/plotToFile.m
new file mode 100644
index 0000000000000000000000000000000000000000..a5fc72f656e0c723daa6968841873aeeae339e64
--- /dev/null
+++ b/lib/storage/@LevelData/plotToFile.m
@@ -0,0 +1,45 @@
+function plotToFile(obj, xVariable, yVariable)
+%%PLOTTOFILE creates plots of various scalar variables (specified in
+%yVariable) with respect %to the variable xVariable and stores it to a file,
+%filename is generated automatically from information in LevelData object
+% PLOTTOFILE(obj, xVariable, yVariable)
+%
+% See also LevelData/plot
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ xVariable {mustBeTextScalar}
+ end
+
+ arguments (Repeating)
+ yVariable {mustBeTextScalar}
+ end
+
+ if isempty(yVariable)
+ yVariable = setdiff(obj.label, xVariable);
+ end
+
+ ensureFolderExists(obj.foldername);
+ h = createStandardFigure();
+ obj.plot(xVariable, yVariable{:});
+
+ % Export plot
+ print(h, '-dpng', '-r600', obj.foldername + '/' + obj.filename + '.png');
+ close(h);
+end
diff --git a/lib/storage/@LevelData/plotTriangulation.m b/lib/storage/@LevelData/plotTriangulation.m
new file mode 100644
index 0000000000000000000000000000000000000000..caae9243368c7b71ea42969b4b75d276acc8a8c0
--- /dev/null
+++ b/lib/storage/@LevelData/plotTriangulation.m
@@ -0,0 +1,48 @@
+function ax = plotTriangulation(obj, jLevel, opt)
+%%PLOTTRIANGULATION creates (a series of) plots for triangulations (must be
+%stored as pair c4n (double: nNodes x 2) and n4e (int: nElem x 3) or as
+%pair coordinates (double: 2 x nNodes) and elements (int: 3 x nElem)
+% ax = PLOTTRIANGULATION(obj) plots all levels
+% ax = PLOTTRIANGULATION(obj, jLevel) plots levels specified in jLevel
+% ax = PLOTTRIANGULATION(obj, jLevel, opt) additionally specifies options:
+% - timePerLevel (default 1s)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jLevel (1,:) {mustBeIndexVector} = ':'
+ opt.timePerLevel {mustBeNumeric} = 1
+ end
+
+ if strcmp(jLevel, ':')
+ jLevel = 1:obj.nLevel;
+ end
+
+ % Create axes object
+ ax = gca;
+
+ % Plot every frame
+ % TODO: store whole mesh such that also boundary can be plotted?
+ for k = jLevel
+ coordinates = obj.level(k).coordinates;
+ elements = obj.level(k).elements;
+ mesh = Mesh(coordinates, elements, {});
+ plot(mesh, 'FaceAlpha', 0.0);
+ pause(opt.timePerLevel);
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/plotTriangulationToFile.m b/lib/storage/@LevelData/plotTriangulationToFile.m
new file mode 100644
index 0000000000000000000000000000000000000000..65db0d8d202941a4609e4ada808382c5bfe55e64
--- /dev/null
+++ b/lib/storage/@LevelData/plotTriangulationToFile.m
@@ -0,0 +1,79 @@
+function plotTriangulationToFile(obj, jLevel, opt)
+%%PLOTTRIANGULATIONTOFILE creates (a series of) plots for triangulations
+%and stores it to a file, filename is generated automatically from
+%information in LevelData object
+% ax = PLOTTRIANGULATIONTOFILE(obj) plots all levels
+% ax = PLOTTRIANGULATIONTOFILE(obj, jLevel) plots levels specified in jLevel
+% ax = PLOTTRIANGULATIONTOFILE(obj, jLevel, opt) additionally specifies options:
+% - timePerLevel (default 1s)
+%
+% See also LevelData/plotTriangulation
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jLevel (1,:) {mustBeIndexVector} = ':'
+ opt.timePerLevel {mustBeNumeric} = 1
+ end
+
+ if strcmp(jLevel, ':')
+ jLevel = 1:obj.nLevel;
+ end
+
+ ensureFolderExists(obj.foldername);
+ h = createStandardFigure();
+
+ if length(jLevel) == 1
+ % Create single image
+ coordinates = obj.level(jLevel).coordinates;
+ elements = obj.level(jLevel).elements;
+ mesh = Mesh(coordinates, elements, {});
+ plot(mesh);
+ print(h, '-dpng', '-r600', obj.foldername + '/' + obj.filename + '.png');
+ else
+ % Create video as series of images
+ createTriangulationAVI(obj, h, jLevel, opt.timePerLevel);
+ end
+
+ close(h);
+end
+
+
+function createTriangulationAVI(obj, h, jLevel, secondsPerFrame)
+ % video quality parameter
+ FPS = 1 / secondsPerFrame;
+
+ % Create video writer
+ videowriter = VideoWriter(obj.foldername + '/' + obj.filename + '.avi');
+ set(videowriter, 'FrameRate', FPS);
+ open(videowriter);
+
+ ax = gca();
+ set(ax, 'NextPlot', 'replaceChildren');
+
+ % Plot every frame
+ for k = jLevel
+ coordinates = obj.level(k).coordinates;
+ elements = obj.level(k).elements;
+ mesh = Mesh(coordinates, elements, {});
+ plot(mesh, 'FaceAlpha', 0.0);
+ writeVideo(videowriter, getframe(h));
+ end
+
+ close(videowriter);
+end
diff --git a/lib/storage/@LevelData/printHeader.m b/lib/storage/@LevelData/printHeader.m
new file mode 100644
index 0000000000000000000000000000000000000000..3888afcbde7a1d02ca10127356b141b45f0e5303
--- /dev/null
+++ b/lib/storage/@LevelData/printHeader.m
@@ -0,0 +1,32 @@
+function printHeader(obj)
+%%PRINTHEADER prints header for output to command line
+% PRINTHEADER(obj)
+%
+% See also LevelData/printLevel
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ % Create header with variable names
+ specifier = obj.getHeaderSpecifier(' ');
+ header = sprintf(specifier, obj.scalarVariable{:});
+
+ % Print header to command line
+ fprintf(header);
+
+ % Print separating horizontal line
+ fprintf([repmat('-', 1, length(header)-1), '\n']);
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/printLevel.m b/lib/storage/@LevelData/printLevel.m
new file mode 100644
index 0000000000000000000000000000000000000000..110e49f1c597600ba3cedd711f4b866ad4497d72
--- /dev/null
+++ b/lib/storage/@LevelData/printLevel.m
@@ -0,0 +1,49 @@
+function printLevel(obj, jLevel)
+%%PRINTLEVEL prints information on the given levels to command line
+% PRINTLEVEL(obj) defaults to final level
+% PRINTLEVEL(obj, jLevel)
+%
+% See also LevelData/printHeader
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jLevel (1,:) double = obj.nLevel
+ end
+
+ % Print header in case of plotting the first level
+ if obj.isInitialRun()
+ obj.printHeader();
+ end
+
+ % Iterate over given list of levels
+ specifier = obj.getFormatSpecifier(' ');
+ for k = 1:length(jLevel)
+ % Extract data of all variables
+ data = cell(obj.nVariable, 1);
+ ind = 0;
+ for jVariable = 1:obj.nVariable
+ if obj.isScalar(jVariable)
+ ind = ind + 1;
+ data{ind} = obj.level(jLevel(k)).(obj.label{jVariable});
+ end
+ end
+ % Print information on current level to command line
+ fprintf(specifier, data{1:ind});
+ end
+end
diff --git a/lib/storage/@LevelData/saveToFile.m b/lib/storage/@LevelData/saveToFile.m
new file mode 100644
index 0000000000000000000000000000000000000000..d18503ed503331bff2e8355fb61604c082f80370
--- /dev/null
+++ b/lib/storage/@LevelData/saveToFile.m
@@ -0,0 +1,35 @@
+function saveToFile(obj, folder, file)
+%%SAVETOFILE saves this object to a file with automatically
+%generated folder and file names if no optional arguments are
+%specified
+% obj.SAVETOFILE()
+% obj.SAVETOFILE(folder, file)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ folder {mustBeTextScalar} = obj.foldername
+ file {mustBeTextScalar} = obj.filename
+ end
+
+ % Create problem- and method-specific folder
+ ensureFolderExists(folder);
+ % Save this object to file
+ pathToMat = folder + '/' + file + '.mat';
+ save(pathToMat, 'obj', '-v7.3');
+end
diff --git a/lib/storage/@LevelData/saveToTable.m b/lib/storage/@LevelData/saveToTable.m
new file mode 100644
index 0000000000000000000000000000000000000000..7527c78ce0d47270bd01b9d2c35d1efedcb504ba
--- /dev/null
+++ b/lib/storage/@LevelData/saveToTable.m
@@ -0,0 +1,48 @@
+function saveToTable(obj, separator)
+%%SAVETOTABLE stores all level information in a text file, defaults to
+%csv-format (separator = ','), filename is generated automatically from
+%information in LevelData object
+% SAVETOTABLE(obj)
+% SAVETOTABLE(obj, separator)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ separator {mustBeTextScalar} = ','
+ end
+
+ % Open file
+ ensureFolderExists(obj.foldername);
+ fid = fopen(obj.foldername + '/' + obj.filename + '.csv', 'w');
+
+ % Save header to file
+ specifier = obj.getHeaderSpecifier(separator);
+ fprintf(fid, specifier, obj.scalarVariable{:});
+
+ % Save information on each level to file
+ specifier = obj.getFormatSpecifier(separator);
+ for k = 1:obj.nLevel
+ data = cell(obj.nScalarVariable, 1);
+ for j = 1:obj.nScalarVariable
+ data{j} = obj.level(k).(obj.scalarVariable{j});
+ end
+ fprintf(fid, specifier, data{:});
+ end
+
+ fclose(fid);
+end
diff --git a/lib/storage/@LevelData/set.m b/lib/storage/@LevelData/set.m
new file mode 100644
index 0000000000000000000000000000000000000000..7a3ee72195dcda906172a451be472995e9f035ab
--- /dev/null
+++ b/lib/storage/@LevelData/set.m
@@ -0,0 +1,65 @@
+function set(obj, jLevel, variableName, value)
+%%SET stores specified data to this LevelData object for the specified list
+%jLevel of levels. The variableName/value pairs can be repeating, their first
+%dimension must correspond to number of levels
+% SET(obj, jLevel, variableName, value, ...)
+% SET(obj, ':', variableName, value, ...)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jLevel {mustBeIndexVector} = ':'
+ end
+
+ arguments (Repeating)
+ variableName {mustBeTextScalar}
+ value
+ end
+
+ if strcmp(jLevel, ':')
+ jLevel = 1:obj.nLevel;
+ end
+
+ % Determine number of arguments
+ nArgument = length(variableName);
+
+ % Store data to LevelData object
+ for j = 1:nArgument
+ name = variableName{j};
+ val = value{j};
+
+ if isa(val, 'ValueDetails')
+ % Copy type argument
+ currentType = val;
+ valueList = repmat({nan}, length(jLevel), 1);
+ else
+ [valueList, currentType] = splitIntoLevelWiseData(length(jLevel), val);
+ end
+
+ % Store type and category
+ if ~ismember(name, obj.label)
+ obj.type(name) = currentType.adaptPrintWidthTo(name);
+ obj.category(name) = DataCategory.ERROR;
+ end
+
+ % Store level-oriented data
+ for k = 1:length(jLevel)
+ obj.level(jLevel(k)).(name) = valueList{k};
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/setAbsolute.m b/lib/storage/@LevelData/setAbsolute.m
new file mode 100644
index 0000000000000000000000000000000000000000..d5fbf92a22bc227a71c503f2f1e56d0cd875aa23
--- /dev/null
+++ b/lib/storage/@LevelData/setAbsolute.m
@@ -0,0 +1,68 @@
+function setAbsolute(obj, jLevel, variableName, value)
+%%SETABSOLUTE stores specified data of absolute-niveau type (constants /
+%eigenvalues / efficiency indices / ...) to this LevelData object for the
+%specified list jLevel of levels, repeated pairs of variable names and values,
+%first dimension of data must correspond to number of levels
+% SETABSOLUTE(obj, jLevel, variableName, value, ...)
+% SETABSOLUTE(obj, ':', variableName, value, ...)
+%
+% See also LevelData/set
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jLevel {mustBeIndexVector} = ':'
+ end
+
+ arguments (Repeating)
+ variableName {mustBeTextScalar}
+ value
+ end
+
+ if strcmp(jLevel, ':')
+ jLevel = 1:obj.nLevel;
+ end
+
+ % Determine number of arguments
+ nArgument = length(variableName);
+
+ % Store data to LevelData object
+ for j = 1:nArgument
+ name = variableName{j};
+ val = value{j};
+
+ if isa(val, 'ValueDetails')
+ % Copy type argument
+ currentType = val;
+ valueList = repmat({nan}, length(jLevel), 1);
+ else
+ [valueList, currentType] = splitIntoLevelWiseData(length(jLevel), val);
+ end
+
+ % Store type and category
+ if ~ismember(name, obj.absoluteVariable)
+ obj.type(name) = currentType.adaptPrintWidthTo(name);
+ obj.category(name) = DataCategory.ABSOLUTE;
+ end
+
+ % Store level-oriented data
+ for k = 1:length(jLevel)
+ obj.level(jLevel(k)).(name) = valueList{k};
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelData/setTime.m b/lib/storage/@LevelData/setTime.m
new file mode 100644
index 0000000000000000000000000000000000000000..c492f28de989cb54686235506c863c0ebceeb1ed
--- /dev/null
+++ b/lib/storage/@LevelData/setTime.m
@@ -0,0 +1,68 @@
+function setTime(obj, jLevel, variableName, value)
+%%SETTIME stores specified data of increasing type (time / ndof / ...) to
+%this LevelData object for the specified list jLevel of levels, repeated pairs
+%of variable names and values, first dimension of data must correspond to number
+%of levels
+% SETTIME(obj, jLevel, variableName, value, ...)
+% SETTIME(obj, ':', variableName, value, ...)
+%
+% See also LevelData/set
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jLevel {mustBeIndexVector} = ':'
+ end
+
+ arguments (Repeating)
+ variableName {mustBeTextScalar}
+ value
+ end
+
+ if strcmp(jLevel, ':')
+ jLevel = 1:obj.nLevel;
+ end
+
+ % Determine number of arguments
+ nArgument = length(variableName);
+
+ % Store data to LevelData object
+ for j = 1:nArgument
+ name = variableName{j};
+ val = value{j};
+
+ if isa(val, 'ValueDetails')
+ % Copy type argument
+ currentType = val;
+ valueList = repmat({nan}, length(jLevel), 1);
+ else
+ [valueList, currentType] = splitIntoLevelWiseData(length(jLevel), val);
+ end
+
+ % Store type and category
+ if ~ismember(name, obj.label)
+ obj.type(name) = currentType.adaptPrintWidthTo(name);
+ obj.category(name) = DataCategory.TIME;
+ end
+
+ % Store level-oriented data
+ for k = 1:length(jLevel)
+ obj.level(jLevel(k)).(name) = valueList{k};
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelDataCollection/LevelDataCollection.m b/lib/storage/@LevelDataCollection/LevelDataCollection.m
new file mode 100644
index 0000000000000000000000000000000000000000..bf886751ea8dee3113f5e910143440ba1514010c
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/LevelDataCollection.m
@@ -0,0 +1,182 @@
+classdef LevelDataCollection < handle
+%%LEVELDATACOLLECTION class representing a collection of LevelData objects
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ properties
+ metaData
+ % Root folder for file storage
+ root = 'results'
+ % Cell array of LevelData objects
+ item = cell(0)
+ end
+
+ properties (Dependent)
+ % Path to folder for file storage
+ foldername (1,:) char
+ % Name of file for storage
+ filename (1,:) char
+ % Number of items
+ nItem (1,1) integer
+ % Cell array of names of timing variables
+ timeVariable (1,:) cell
+ % Number of time variables
+ nTimeVariable (1,1) int32
+ % Boolean indicating whether data has been stored for a single level
+ isInitialRun (1,1) boolean
+ end
+
+ properties (Access=protected)
+ % Separator for printing to command line
+ separator
+ % Minimal width of variables for printing in command line
+ minimalWidth = 8;
+ end
+
+ methods
+ %% CONSTRUCTOR
+ function obj = LevelDataCollection(rootpath, identifier)
+ %%LEVELDATACOLLECTION
+
+ % Set root path for file storage
+ if nargin >= 1
+ obj.root = rootpath;
+ end
+
+ obj.metaData = dictionary(...
+ "problem", "problem", ...
+ "domain", "domain", ...
+ "method", "method", ...
+ "identifier", "main", ...
+ "hostname", string(getHostname()), ...
+ "timestamp", string(datetime('now', 'Format', 'yyyy-MM-dd_HH:mm:ss')));
+
+ % Set identifier
+ if nargin >= 2
+ obj.metaData("identifier") = identifier;
+ end
+ end
+
+ %% GET FUNCTIONS FOR DEPENDENT VARIABLES
+ function nItem = get.nItem(obj)
+ nItem = length(obj.item);
+ end
+
+ function timeVariable = get.timeVariable(obj)
+ if isempty(obj.item)
+ timeVariable = {};
+ else
+ timeVariable = obj.item{1}.timeVariable;
+ end
+ end
+
+ function nVariable = get.nTimeVariable(obj)
+ nVariable = length(obj.timeVariable);
+ end
+
+ function bool = get.isInitialRun(obj)
+ bool = (obj.nItem <= 1);
+ end
+
+ function path = get.foldername(obj)
+ rootpath = obj.root;
+ if numel(rootpath) > 0 && ~strcmp(rootpath(end), '/')
+ rootpath = [rootpath, '/'];
+ end
+ path = rootpath + strjoin([obj.metaData("problem"), obj.metaData("domain"), obj.metaData("method")], '_');
+ end
+
+ function file = get.filename(obj)
+ file = obj.metaData("identifier") + '_collection';
+ end
+
+ function spec = getHeaderSpecifier(obj, separator)
+ % Creates formatting string for the header of the output to command line
+ spec = cell(1, obj.nTimeVariable+1);
+ spec{1} = ' run';
+ for j = 1:obj.nTimeVariable
+ t = obj.item{1}.type(obj.timeVariable{j});
+ spec{j+1} = assembleSpecifier(obj.getWidth(t), 's');
+ end
+ spec = strjoin(spec, separator) + "\n";
+ end
+
+ function spec = getFormatSpecifier(obj, separator)
+ % Creates formatting string for printing to command line
+ spec = cell(1, obj.nTimeVariable+1);
+ spec{1} = '%4d';
+ for j = 1:obj.nTimeVariable
+ t = obj.item{1}.type(obj.timeVariable{j});
+ spec{j+1} = assembleSpecifier(obj.getWidth(t), t.formatSpec);
+ end
+ spec = strjoin(spec, separator) + "\n";
+ end
+
+ %% READ ITEM DATA
+ output = get(obj, jItem, varargin)
+
+ %% MODIFY ITEMS
+ set(obj, jItem, varargin)
+
+ function append(obj, data)
+ %%APPEND simplified addition specified of data to this list, one or
+ %more LevelData objects
+ % APPEND(obj, data)
+
+ arguments
+ obj
+ end
+
+ arguments (Repeating)
+ data LevelData
+ end
+
+ indices = obj.nItem + (1:length(data));
+ obj.set(indices, data{:});
+ end
+
+ function remove(obj, indices)
+ %%REMOVE removes the specified list of objects
+ % REMOVE(obj, indices)
+ indices = intersect(indices, 1:obj.nItem);
+ obj.item{indices} = [];
+ end
+
+ %% OUTPUT ITEM DATA
+ printHeader(obj)
+ printItem(obj, jItem)
+ printStatistics(obj, variableName)
+ plotStatistics(obj, xVariable, variableName)
+
+ %% EXPORT COMPLETE DATA
+ saveToFile(obj, folder, file)
+ saveToTable(obj, separator)
+ end
+
+ methods (Access = private)
+ function width = getWidth(obj, type)
+ width = num2str(max(type.printWidth, obj.minimalWidth));
+ end
+
+ printTable(obj, fid, variableName, data, separator)
+ end
+end
+
+function spec = assembleSpecifier(width, format)
+ spec = ['%', num2str(width), format];
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelDataCollection/get.m b/lib/storage/@LevelDataCollection/get.m
new file mode 100644
index 0000000000000000000000000000000000000000..025502e872f02d5df631df4b34724f74325186e6
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/get.m
@@ -0,0 +1,60 @@
+function output = get(obj, jItem, variableName)
+%%GET extracts data from this LevelDataCollection object for a given list
+%jItem of item numbers. One or more variables to be returned can be specified by
+%their names.
+% output = GET(obj, jItem, variableName, ...)
+% output = GET(obj, ':', variableName, ...)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jItem {mustBeIndexVector} = ':'
+ end
+
+ arguments (Repeating)
+ variableName {mustBeTextScalar}
+ end
+
+ if strcmp(jItem, ':')
+ jItem = 1:obj.nItem;
+ end
+
+ % Initialise output variable
+ output = cell(1, length(variableName));
+
+ % Determine number of levels (uses first item!)
+ nLevel = obj.item{1}.nLevel;
+
+ % Iterate over variables
+ for jVariable = 1:length(variableName)
+ name = variableName{jVariable};
+ % Check for presence of variable
+ if ~ismember(name, obj.timeVariable)
+ warning(['Variable ', name, ' not found']);
+ data = [];
+ else
+ data = zeros(nLevel, length(jItem));
+ % Iterate over items
+ for k = 1:length(jItem)
+ data(:,k) = obj.item{jItem(k)}.get(':', name);
+ end
+ end
+ % Store to return variable
+ output{jVariable} = data;
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelDataCollection/plotStatistics.m b/lib/storage/@LevelDataCollection/plotStatistics.m
new file mode 100644
index 0000000000000000000000000000000000000000..3c81e352890e4e5efae7738d45e13f209b75d13c
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/plotStatistics.m
@@ -0,0 +1,90 @@
+function ax = plotStatistics(obj, xVariable, yVariable)
+%%PLOTSTATISTICS plots with statistical information (mean value with min
+%and max value) of scalar time variables (specified in yVariable) with
+%respect to the variable xVariable, returns handle to axis object
+% PLOTSTATISTICS(obj, xVariable, yVariable, ...)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ xVariable {mustBeTextScalar}
+ end
+
+ arguments (Repeating)
+ yVariable
+ end
+
+ % Choose time variables for plotting
+ if isempty(yVariable)
+ yVariable = obj.timeVariable;
+ else
+ yVariable = intersect(yVariable, obj.timeVariable);
+ end
+
+ % Extract data for x-axis
+ xData = obj.item{1}.get(':', xVariable);
+
+ % Extract data from each item and on each level
+ data = obj.get(':', yVariable{:});
+
+ % Create plot
+ ax = plotData(xData, yVariable, data);
+
+ % Add title
+ title(ax, 'Runtime plot');
+
+ % Add axes labels
+ xlabel(ax, xVariable);
+ ylabel(ax, 'runtime (in sec.)');
+
+ % Update legend
+ configureLegend(ax, 'northwest');
+end
+
+
+function ax = plotData(xData, variableName, data)
+%%PLOTDATA auxiliary private function for creation of plots, returns handle
+%to axis object
+% ax = PLOTDATA(xData, variableName, data)
+
+ % Create handle to currently active axis object
+ ax = gca;
+
+ % Create unified list for formatting plot lines
+ [COLOURORDER, MARKER] = getPlotStyle();
+
+ % Iterate over given variables
+ for j = 1:length(variableName)
+ % Compute statistics
+ yMean = mean(data{j}, 2);
+ yMin = min(data{j}, [], 2);
+ yMax = max(data{j}, [], 2);
+ % Extract label for legend from dictionary
+ variableLabel = variableName{j};
+ % Create plot
+ errorbar(ax, xData, yMean, ...
+ yMean - yMin, yMax - yMean, ...
+ 'LineStyle', '-', 'Marker', '.', ...
+ 'Color', COLOURORDER(j,:), ...
+ 'DisplayName', variableLabel);
+ % Use double-logarithmic scales
+ set(ax, 'XScale', 'log', 'YScale', 'log');
+ % Add new line into the current figure when calling plotConvergence again
+ hold(ax, 'on');
+ end
+end
diff --git a/lib/storage/@LevelDataCollection/printHeader.m b/lib/storage/@LevelDataCollection/printHeader.m
new file mode 100644
index 0000000000000000000000000000000000000000..b3225a5b71aa7b68a27a5c8de60bdfc82d0324d9
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/printHeader.m
@@ -0,0 +1,32 @@
+function printHeader(obj)
+%%PRINTHEADER prints header for output to command line
+% PRINTHEADER(obj)
+%
+% See also LevelDataCollection/printItem
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ % Create header with variable names
+ header = sprintf(obj.getHeaderSpecifier(' '), obj.timeVariable{:});
+
+ % Print header to command line
+ fprintf(header);
+
+ % Print separating horizontal line
+ fprintf([repmat('-', 1, length(header)-1), '\n']);
+end
diff --git a/lib/storage/@LevelDataCollection/printItem.m b/lib/storage/@LevelDataCollection/printItem.m
new file mode 100644
index 0000000000000000000000000000000000000000..b83bfd210266fae755f8f4acf813b735dafee24e
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/printItem.m
@@ -0,0 +1,48 @@
+function printItem(obj, jItem)
+%%PRINTITEM prints information on the given list of item indices to command
+%line
+% PRINTITEM(obj) defaults to final level
+% PRINTITEM(obj, jItem)
+%
+% See also LevelDataCollection/printHeader
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jItem = obj.nItem
+ end
+
+ % Print header in case of plotting the first level
+ if obj.isInitialRun()
+ obj.printHeader();
+ end
+
+ % Iterate over given list of item indices
+ for k = 1:length(jItem)
+ % Extract data of all time variables
+ data = cell(obj.nTimeVariable+1, 1);
+ data{1} = jItem(k);
+ ind = 1;
+ for l = 1:obj.nTimeVariable
+ ind = ind + 1;
+ data{ind} = obj.item{jItem(k)}.level.(obj.timeVariable{l});
+ end
+ % Print information on current item to command line
+ fprintf(obj.getFormatSpecifier(' '), data{1:ind});
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelDataCollection/printStatistics.m b/lib/storage/@LevelDataCollection/printStatistics.m
new file mode 100644
index 0000000000000000000000000000000000000000..1c9766a445b5d2e4ddf419950995a9145ce54850
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/printStatistics.m
@@ -0,0 +1,50 @@
+function printStatistics(obj, variable)
+%%PRINTSTATISTICS prints statististical information of scalar time
+%variables (specified in variable)
+% PRINTSTATISTICS(obj, variable, ...)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ end
+
+ arguments (Repeating)
+ variable {mustBeTextScalar}
+ end
+
+ % Choose time variables for plotting
+ if isempty(variable)
+ variable = obj.timeVariable;
+ else
+ variable = intersect(variable, obj.timeVariable);
+ end
+
+ % Extract data from each item and on each level
+ data = obj.get(':', variable{:});
+
+ % Print separator between tables
+ fprintf('\n');
+
+ % Print table for each given variable
+ for j = 1:length(variable)
+ obj.printTable(1, variable{j}, data{j}, ' ');
+ fprintf('\n\n');
+ end
+end
+
+
diff --git a/lib/storage/@LevelDataCollection/printTable.m b/lib/storage/@LevelDataCollection/printTable.m
new file mode 100644
index 0000000000000000000000000000000000000000..7756663e2784ed9bde0f747d0d9a34c3d595f3f4
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/printTable.m
@@ -0,0 +1,48 @@
+function printTable(obj, fid, variableName, data, separator)
+%%PRINTTABLE auxiliary private function for printing statististical
+%information for data of scalar time variables variableName to command
+%line (fid=1) or file (specified by file identifier fid)
+% PRINTTABLE(obj, fid, variableName, data)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ % Define formatting strings for title and headline
+ TITLE = ['# TIME STATISTICS - ', variableName];
+ HEADLINE = sprintf(['%5s', repmat([separator, '%11s'], 1, 5)],...
+ 'level', 'mean', 'median', 'std', 'min', 'max');
+
+ if fid == 1
+ % Print to command line (including title)
+ fprintf(fid, [TITLE, '\n\n', HEADLINE, '\n', ...
+ repmat('-', 1, length(HEADLINE)), '\n']);
+ else
+ % Print to file
+ fprintf(fid, [HEADLINE, '\n']);
+ end
+
+ % Compute statistical information
+ statisticalData = [mean(data, 2),...
+ median(data, 2),...
+ std(data, 0, 2),...
+ min(data, [], 2),...
+ max(data, [], 2)];
+
+ % Print information
+ fprintf(fid, ['%5d', repmat([separator, '%8.5e'], 1, 5), '\n'],...
+ [1:size(data, 1); permute(statisticalData, [2 1])]);
+end
\ No newline at end of file
diff --git a/lib/storage/@LevelDataCollection/saveToFile.m b/lib/storage/@LevelDataCollection/saveToFile.m
new file mode 100644
index 0000000000000000000000000000000000000000..d1aca0c2a977d9513ad522b6b3d40cb21047c62b
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/saveToFile.m
@@ -0,0 +1,35 @@
+function saveToFile(obj, folder, file)
+%%SAVETOFILE saves this object to a file with automatically
+%generated folder and file names if no optional arguments are
+%specified
+% obj.SAVETOFILE()
+% obj.SAVETOFILE(folder, file)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ folder = obj.foldername
+ file = obj.filename
+ end
+
+ % Create problem- and method-specific folder
+ ensureFolderExists(folder);
+
+ % Save this object to file
+ save(folder + '/' + file + '.mat', 'obj', '-v7.3');
+end
diff --git a/lib/storage/@LevelDataCollection/saveToTable.m b/lib/storage/@LevelDataCollection/saveToTable.m
new file mode 100644
index 0000000000000000000000000000000000000000..1ec56e3af4dd4af377b38c19933271e5e15dd937
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/saveToTable.m
@@ -0,0 +1,40 @@
+function saveToTable(obj, separator)
+%%SAVETOTABLE stores statistical information for all time variables to a
+%text file, defaults to csv-format (separator = ','), filename is generated
+%automatically from information in LevelData object
+% SAVETOTABLE(obj)
+% SAVETOTABLE(obj, separator)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ separator {mustBeTextScalar} = ','
+ end
+
+ % Create problem- and method-specific folder
+ ensureFolderExists(obj.foldername);
+
+ % Save data for each variable to a separate file
+ data = obj.get(':', obj.timeVariable{:});
+ for j = 1:obj.nTimeVariable
+ fid = fopen(obj.foldername + '/' + obj.filename + '_' + ...
+ obj.timeVariable{j} + '.csv', 'w');
+ obj.printTable(fid, obj.timeVariable{j}, data{j}, separator);
+ fclose(fid);
+ end
+end
diff --git a/lib/storage/@LevelDataCollection/set.m b/lib/storage/@LevelDataCollection/set.m
new file mode 100644
index 0000000000000000000000000000000000000000..6ca051f4beade0d6e069ee543d982243f2498f45
--- /dev/null
+++ b/lib/storage/@LevelDataCollection/set.m
@@ -0,0 +1,39 @@
+function set(obj, jItem, data)
+%%SET stores specified list of LevelData objects to the indices determined
+%in the jItem array
+% SET(obj, jItem, data)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ obj
+ jItem
+ end
+
+ arguments (Repeating)
+ data LevelData
+ end
+
+ % Number of items and item indices must coincide
+ assert(length(jItem) == length(data), ...
+ 'Number of indices must equal number of given arguments');
+
+ % Store items
+ for k = 1:length(jItem)
+ obj.item{jItem(k)} = data{k};
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/DataCategory.m b/lib/storage/DataCategory.m
new file mode 100644
index 0000000000000000000000000000000000000000..f4ea739eee9d02b7de5f64161fdc34135701490f
--- /dev/null
+++ b/lib/storage/DataCategory.m
@@ -0,0 +1,23 @@
+classdef DataCategory
+ enumeration
+ ERROR(@loglog, 'Convergence history plot', 'error', 'northeast')
+ ABSOLUTE(@semilogx, 'Value plot', 'value', 'northeast')
+ TIME(@loglog, 'Time plot', 'runtime', 'northwest')
+ end
+
+ properties (SetAccess=immutable)
+ plotFunction function_handle
+ title string
+ yLabel string
+ legendLocation string
+ end
+
+ methods
+ function obj = DataCategory(plotFunction, title, yLabel, legendLocation)
+ obj.plotFunction = plotFunction;
+ obj.title = title;
+ obj.yLabel = yLabel;
+ obj.legendLocation = legendLocation;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/TimeIt.m b/lib/storage/TimeIt.m
new file mode 100644
index 0000000000000000000000000000000000000000..dfc7b06e12aff8d96cb7c38535242a70f3378a3d
--- /dev/null
+++ b/lib/storage/TimeIt.m
@@ -0,0 +1,77 @@
+function leveldatacollection = TimeIt(identifier, nRun, storeResults, functionName, arguments)
+%%TIMEIT function wrapper for multiple runs of functions storing timing
+%data in LevelData objects
+% leveldatacollection = TIMEIT(identifier, nRun, functionName, arguments, ...)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ arguments
+ identifier {mustBeTextScalar}
+ nRun (1,1) double
+ storeResults (1,1) logical
+ functionName {mustBeTextScalar}
+ end
+
+ arguments (Repeating)
+ arguments
+ end
+
+ % Welcome statement
+ fprintf('\n#\n# MooAFEM - TIME MEASUREMENT\n');
+ fprintf('# Current Time: %s\n#\n\n', datestr(now));
+ fprintf('This is TimeIt wrapper for function "%s"\n\n', functionName);
+
+ % Initialisation of collection for LevelData objects
+ leveldatacollection = LevelDataCollection();
+ leveldatacollection.metaData("identifier") = identifier;
+
+ % Run experiments
+ for j = 1:nRun
+ % Run current experiment
+ temporaryIdentifier = sprintf('timingRun%04d', j);
+ outputList = fevalc(functionName, arguments);
+ leveldata = outputList{1};
+ leveldata.metaData("identifier") = temporaryIdentifier;
+ % Remove unused information
+ leveldata.removeNonscalar();
+ % Store information
+ leveldatacollection.append(leveldata);
+ if storeResults
+ % Save intermediate collection to file
+ leveldatacollection.saveToFile();
+ end
+ % Print information on current run
+ leveldatacollection.printItem();
+ end
+
+ % Closing
+ fprintf('\n');
+end
+
+
+function output = fevalc(functionName, arguments) %#ok<INUSD>
+%%FEVALC suppresses output to commandline
+
+ % Create function call
+ functioncall = 'feval(functionName, arguments{:})';
+ % Call function
+ [~, output] = evalc(functioncall);
+ % Store output in cell variable
+ if ~iscell(output)
+ output = {output};
+ end
+end
diff --git a/lib/storage/Type/RawType.m b/lib/storage/Type/RawType.m
new file mode 100644
index 0000000000000000000000000000000000000000..94bdffa19c127733384cf6dbebfe85e9a00fd95f
--- /dev/null
+++ b/lib/storage/Type/RawType.m
@@ -0,0 +1,18 @@
+classdef RawType
+ enumeration
+ INT('d')
+ FLOAT('.5e')
+ TEXT('s')
+ end
+
+ properties (SetAccess=immutable, GetAccess=public)
+ formatSpec
+
+ end
+
+ methods
+ function obj = RawType(formatSpec)
+ obj.formatSpec = formatSpec;
+ end
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/Type/ValueDetails.m b/lib/storage/Type/ValueDetails.m
new file mode 100644
index 0000000000000000000000000000000000000000..e82680e5e5bc514d1530283422e44ee2fadaa34d
--- /dev/null
+++ b/lib/storage/Type/ValueDetails.m
@@ -0,0 +1,107 @@
+classdef ValueDetails
+%%ValueDetails class representing the type and printing parameters of a level
+%oriented data
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ properties (GetAccess=public, SetAccess=protected)
+ % Raw type
+ rawType
+ % Shape of data (scalar, vector, or array)
+ shape
+ % Number of characters required for printing
+ printWidth
+ end
+
+ properties (Dependent)
+ % Format specification for printing
+ formatSpec
+ % Boolean determining if type is scalar
+ isScalar (1,1) boolean
+ end
+
+ methods
+ %% CONSTRUCTOR
+ function obj = ValueDetails(rawType, shape, printWidth)
+ %%VALUEDETAILS creates an instance of ValueDetails given an example
+ %value of the corresponding data
+ % obj = VALUEDETAILS(rawType, shape, printWidth)
+
+ obj.rawType = rawType;
+ obj.shape = shape;
+ obj.printWidth = printWidth;
+ end
+
+ function obj = adaptPrintWidthTo(obj, name)
+ %ADAPTPRINTWIDTHTO updates the print width considering the length of
+ %variable name
+ obj.printWidth = max(obj.printWidth, length(name));
+ end
+
+ function formatSpec = get.formatSpec(obj)
+ formatSpec = obj.rawType.formatSpec;
+ end
+
+ function bool = get.isScalar(obj)
+ bool = strcmp(obj.shape, 's');
+ end
+ end
+
+ % AUXILIARY
+ methods (Static)
+ function details = fromExample(value)
+ %%FROMEXAMPLE creates an instance of ValueDetails given an example
+ %value of the corresponding data
+ % obj = FROMEXAMPLE(value)
+
+ % Determine data type and width
+ if isnan(value)
+ % NaN is treated as integer to avoid problems with plotting
+ rawType = RawType.INT;
+ printWidth = 3;
+ elseif isinteger(value)
+ rawType = RawType.INT;
+ S = cellstr(num2str(value));
+ printWidth = max(cellfun(@(x) length(x), S(:)));
+ elseif isfloat(value)
+ rawType = RawType.FLOAT;
+ floatPrecision = str2double(RawType.FLOAT.formatSpec(2));
+ printWidth = 6 + floatPrecision;
+ elseif ischar(value)
+ rawType = RawType.TEXT;
+ printWidth = length(value);
+ elseif isstring(value)
+ rawType = RawType.TEXT;
+ printWidth = max(cellfun(@(x) length(x), value(:)));
+ else
+ error('Invalid data type');
+ end
+
+ % Determine shape
+ if isscalar(value)
+ shape = 's';
+ elseif isvector(value)
+ shape = 'v';
+ else
+ shape = 'a';
+ end
+
+ details = ValueDetails(rawType, shape, printWidth);
+ end
+ end
+end
diff --git a/lib/storage/Type/splitIntoLevelWiseData.m b/lib/storage/Type/splitIntoLevelWiseData.m
new file mode 100644
index 0000000000000000000000000000000000000000..d399959292ad181a627cda0b7569d19a40c02ef8
--- /dev/null
+++ b/lib/storage/Type/splitIntoLevelWiseData.m
@@ -0,0 +1,64 @@
+
+function [valueList, type] = splitIntoLevelWiseData(nLevel, value)
+%%DETERMINETYPEVALUE auxiliary function for determining type from an array
+%of values
+% [type, valueList] = DETERMINETYPEVALUE(nLevel, variableName, value)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ if nLevel == 1
+ % Array value represents data on a single level
+ type = ValueDetails.fromExample(value);
+ valueList = {value};
+ else
+ if isvector(value)
+ % Array value represents a vector of scalar values for each % level
+ assert(length(value) == nLevel, 'Length of argument invalid');
+ type = ValueDetails.fromExample(value(1));
+ valueList = mat2cell(value(:), ones(nLevel, 1));
+ else
+ % Array value represents an higher dimensional objects for each % level
+ dims = size(value);
+ [levelDim, idx] = extractSingleLevelDimension(dims, nLevel);
+ % Extract values
+ valueList = cell(nLevel, 1);
+ for k = 1:nLevel
+ idx(levelDim) = k;
+ valueList{k} = subsref(value, idx);
+ end
+ type = ValueDetails.fromExample(valueList{end});
+ end
+ end
+end
+
+function [levelDim, idx] = extractSingleLevelDimension(dims, nLevel)
+ levelDim = find(dims == nLevel);
+ if isempty(levelDim)
+ error('Unable to extract data from argument');
+ end
+
+ if length(levelDim) > 1
+ % Several dimensions could represent the level information
+ levelDim = levelDim(1);
+ warning('Unclear dimensions of argument. Assume level-oriented data in dimension %d.', levelDim);
+ end
+
+ % index that is ':' for all dimensions, and where the index in the level
+ % dimension can be changed
+ idx.subs = repmat({':'}, 1, ndims(value));
+ idx.type = '()';
+end
\ No newline at end of file
diff --git a/lib/storage/configureLegend.m b/lib/storage/configureLegend.m
new file mode 100644
index 0000000000000000000000000000000000000000..fe1a7c95992effeec2425a89aab902fb2315572b
--- /dev/null
+++ b/lib/storage/configureLegend.m
@@ -0,0 +1,31 @@
+function configureLegend(ax, location)
+%%CONFIGURELEGEND formats the legend of the current axis and displays it in
+%the specified location
+% CONFIGURELEGEND(ax, location)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ % Create handle to legend in current axis
+ leg = legend(ax);
+
+ % Specify location
+ set(leg, 'Location', location);
+
+ % Set interpreter to latex
+ set(leg, 'Interpreter', 'latex');
+end
diff --git a/lib/storage/createStandardFigure.m b/lib/storage/createStandardFigure.m
new file mode 100644
index 0000000000000000000000000000000000000000..2f1d4cfb83fe436d9b8c68b09ffe1ae40374d046
--- /dev/null
+++ b/lib/storage/createStandardFigure.m
@@ -0,0 +1,26 @@
+function h = createStandardFigure()
+%%CREATESTANDARDFIGURE creates a new figure object with default choice of
+%size, the modification of this file allows a unified formatting of all
+%figures for the output in the LevelData class
+% h = CREATESTANDARDFIGURE()
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ % Create figure
+ h = figure('Units', 'centimeters', 'Position', 3 + [0 0 32 24]);
+
+end
\ No newline at end of file
diff --git a/lib/storage/ensureFolderExists.m b/lib/storage/ensureFolderExists.m
new file mode 100644
index 0000000000000000000000000000000000000000..d29bd1942e064b5e3107b429554e64c9b533e411
--- /dev/null
+++ b/lib/storage/ensureFolderExists.m
@@ -0,0 +1,30 @@
+function ensureFolderExists(path)
+%%ENSUREFOLDEREXISTS creates a folder if it does not exist yet
+% ENSUREFOLDEREXISTS(PATH)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+ if isfolder(path)
+ folderExists = true;
+ else
+ folderExists = mkdir(path);
+ end
+
+ if ~folderExists
+ error('Could not create folder %s', path);
+ end
+end
\ No newline at end of file
diff --git a/lib/storage/getDefaultDictionary.m b/lib/storage/getDefaultDictionary.m
new file mode 100644
index 0000000000000000000000000000000000000000..61e7218ea3474c8c7fdefc169f9d03edb7cc3745
--- /dev/null
+++ b/lib/storage/getDefaultDictionary.m
@@ -0,0 +1,52 @@
+function dict = getDefaultDictionary()
+%%GETDEFAULTDICTIONARY allows to specify default structure arrays for
+%dictionaries for legend labels in plots
+% dict = GETDEFAULTDICTIONARY()
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ % Dictionary for the latex description of data in the
+ % convergence history plot
+ dict =...
+ dictionary('eta' , '$\eta_\ell$',...
+ 'etaVol' , '$\eta_{\textrm{vol}}$',...
+ 'etaJump' , '$\eta_{\textrm{jump}}$',...
+ 'oscDb' , '$\textrm{osc}(\textrm{D} u_{\textrm{D}}, \mathcal{F}_\ell(\Gamma_{\textrm{D}})$',...
+ 'oscNb' , '$\textrm{osc}(t_{\textrm{N}}, \mathcal{F}_\ell(\Gamma_{\textrm{N}})$',...
+ 'oscF' , '$\textrm{osc}(f, \mathcal{T}_\ell)$',...
+ 'mu' , '$|| f - \Pi_0 f ||_{L^2(\Omega)}$',...
+ 'est' , 'built-in',...
+ 'res' , '$LS(f; \sigma_\ell, u_\ell)$',...
+ 'resRed' , '$LS(\Pi_0 f; \sigma_\ell, u_\ell)$',...
+ 'resDiv' , '$|| f + \textrm{div} \sigma_\ell ||_{L^2(\Omega)}$',...
+ 'resL2' , '$|| \sigma_\ell - \nabla u_\ell ||_{L^2(\Omega)}$',...
+ 'resDev' , '$|| \textrm{dev} \sigma_\ell - \textrm{D} u_\ell ||_{L^2(\Omega)}$',...
+ 'resMat' , '$|| C^{-1} \sigma_\ell - \varepsilon u_\ell ||_{L^2(\Omega)}$',...
+ 'errU' , '$|| u - u_\ell ||_{H^1(\Omega)}$',...
+ 'errL2U' , '$|| u - u_\ell ||_{L^2(\Omega)}$',...
+ 'errGradU' , '$||| u - u_\ell |||$',...
+ 'errSigma' , '$|| \sigma - \sigma_\ell ||_{H(\textrm{div},\Omega)}$',...
+ 'errL2Sigma' , '$|| \sigma - \sigma_\ell ||_{L^2(\Omega)}$',...
+ 'errDivSigma' , '$|| \textrm{div} (\sigma - \sigma_\ell) ||_{L^2(\Omega)}$',...
+ 'runtime' , 'runtime',...
+ 'setupTime' , 'setuptime',...
+ 'solutionTime' , 'solutiontime',...
+ 'estimatorTime' , 'estimatortime',...
+ 'refinementTime', 'refinementtime',...
+ 'condition' , 'condition');
+end
diff --git a/lib/storage/getHostname.m b/lib/storage/getHostname.m
new file mode 100644
index 0000000000000000000000000000000000000000..dd77b972d89fc99b6ab719448422de76d5854631
--- /dev/null
+++ b/lib/storage/getHostname.m
@@ -0,0 +1,38 @@
+function hostname = getHostname()
+%%GETHOSTNAME returns the hostname on UNIX systems
+% hostname = GETHOSTNAME()
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ %% READ HOSTNAME
+ if ispc()
+ hostname = 'win';
+ elseif ismac()
+ hostname = 'mac';
+ elseif isunix()
+ [~, hostname] = system('cat /etc/hostname');
+ hostname = hostname(1:end-1);
+ else
+ hostname = 'unknown';
+ end
+
+ %% CHECK RESULT
+ if isempty(hostname)
+ hostname = 'no_hostname';
+ end
+end
diff --git a/lib/storage/getPlotStyle.m b/lib/storage/getPlotStyle.m
new file mode 100644
index 0000000000000000000000000000000000000000..46e0482e8e2e05ac0cee5f42d9cab397e829b3be
--- /dev/null
+++ b/lib/storage/getPlotStyle.m
@@ -0,0 +1,52 @@
+function [colorOrder, marker] = getPlotStyle()
+%%GETPLOTSTYLE creates lists of default colours and markers for plotting
+% [colorOrder, marker] = GETPLOTSTYLE()
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ % Generate a sequence of colours for the graphs
+ % Matlab default
+ colorOrder = [0 0.4470 0.7410;
+ 0.8500 0.3250 0.0980;
+ 0.9290 0.6940 0.1250;
+ 0.4940 0.1840 0.5560;
+ 0.4660 0.6740 0.1880;
+ 0.3010 0.7450 0.9330;
+ 0.6350 0.0780 0.1840];
+
+ % ALTERNATIVE: TU corporate design colours
+ % colorOrder = [0 0.4 0.6;
+ % 0.3922 0.3882 0.3882;
+ % 0 0.4941 0.4431;
+ % 0.7294 0.2745 0.5098;
+ % 0.8824 0.5373 0.1333];
+
+ % ALTERNATIVE: variable length colour sequence
+ % nColor = max(ceil(nVariable^(1/3)), 2);
+ nColor = 3;
+ [red, green, blue] = meshgrid(linspace(0, 1, nColor));
+ % colorOrder = [red(:), green(:), blue(:)];
+ % Append
+ colorOrder = [colorOrder; red(:), green(:), blue(:)];
+ % Remove white as line colour
+ colorOrder = colorOrder(1:end-1,:);
+
+ % Define sequence of markers
+ marker = {'o', '+', 'x', 'square', 'diamond', '*', '^', 'v', ...
+ 'pentagram', 'hexagram', '<', '>', '.', '_', '|'};
+end
diff --git a/runAdditiveSchwarzPreconditioner.m b/runAdditiveSchwarzPreconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..971ffd8d1f8979c73c28c1216fafa76b00187343
--- /dev/null
+++ b/runAdditiveSchwarzPreconditioner.m
@@ -0,0 +1,124 @@
+function leveldata = runAdditiveSchwarzPreconditioner(pMax)
+
+ % Proceed input
+ if nargin < 1
+ pMax = 1e1;
+ end
+
+ % Number of refinement levels
+ nLevels = 2;
+
+ % Initialize LevelData
+ leveldata = LevelData('results');
+ leveldata.problem = 'Poisson';
+ leveldata.domain = 'UnitSquare';
+ leveldata.method = 'Sp_PCG_AdditiveSchwarz';
+
+ % Run experiment for various p
+ for p = 1:pMax
+ % Load mesh
+ mesh = Mesh.loadFromGeometry('unitsquare');
+ %mesh.refineUniform();
+
+ % Create FE space
+ fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', ':');
+ % u = FeFunction(fes);
+
+ % Setup variational problem
+ blf = BilinearForm();
+ lf = LinearForm();
+ blf.a = Constant(mesh, 1);
+ blf.qra = QuadratureRule.ofOrder(max(2*(p-1), 1));
+ lf.f = Constant(mesh, 1);
+
+ % Choose iterative solver
+ solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+ % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
+
+ for k = 1:nLevels-1
+ % Assemble matrices on intermediate meshes
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
+ freeDofs = getFreeDofs(fes);
+ solver.setupSystemMatrix(A(freeDofs,freeDofs));
+ solver.setupRhs(F(freeDofs)); % initializes x0 as well
+ % Mesh refinement
+ mesh.refineUniform();
+ end
+
+ % Assemble matrices on fine mesh
+ A = assemble(blf, fes);
+ if ~issymmetric(A)
+ %warning(['symmetrizing, norm(asym(A)) = ', num2str(norm(A - A', 1))])
+ A = (A + A') / 2;
+ end
+ F = assemble(lf, fes);
+ freeDofs = getFreeDofs(fes);
+ solver.setupSystemMatrix(A(freeDofs,freeDofs));
+ solver.setupRhs(F(freeDofs)); % initializes x0 as well
+
+ % A = A(freeDofs,freeDofs);
+
+ % figure(2);
+ % spy(A);
+ % title('Matrix A');
+
+ % Compute C
+ % C = zeros(length(freeDofs));
+ % for j = 1:length(freeDofs)
+ % ej = zeros(length(freeDofs), 1);
+ % ej(j) = 1;
+ % C(:,j) = solver.preconditionAction(ej);
+ % end
+ % if ~issymmetric(C)
+ % % warning(['symmetrizing, norm(asym(C)) = ', num2str(norm(C - C', 1))])
+ % C = (C + C') / 2;
+ % end
+
+ % figure(3);
+ % CC = C;
+ % CC(abs(C) < 1e-10) = 0;
+ % spy(CC);
+ % title('Matrix C');
+ %
+ % asymC = (C - C');
+
+ % figure(4);
+ % largeAsymC = asymC;
+ % largeAsymC(abs(largeAsymC) < max(abs(largeAsymC(:))) / 10) = 0;
+ % spy(largeAsymC);
+ % title('Large entries of asym(C)');
+
+ % figure(5);
+ % asymC(abs(asymC) < 1e-10) = 0;
+ % spy(asymC);
+ % title('Nonzero entries of asym(C)');
+
+ % Compute conditional number
+ % CA = C*A;
+ % lambda = eig(full(CA));
+ % condition = max(abs(lambda)) / min(abs(lambda));
+ % lambdaMin = eigs(@(x)applyCA(x, A, solver), length(freeDofs), 1, 'smallestreal');
+ % lambdaMax = eigs(@(x)applyCA(x, A, solver), length(freeDofs), 1, 'largestreal');
+ % lambdaMin = eigs(CA, 1, 'smallestabs');
+ % lambdaMax = eigs(CA, 1, 'largestabs');
+ % condition = lambdaMax / lambdaMin;
+
+ % Print result to commandline
+ % leveldata.append('p', uint32(p), 'condition', condition);
+ % leveldata.printLevel();
+ end
+
+ % Plot condition
+ figure();
+ leveldata.plot('p');
+end
+
+
+function CAx = applyCA(x, A, solver)
+ % CAx = solver.preconditionAction(A * x);
+ % CAx = (A * x);
+ CAx = solver.preconditionAction(x);
+end
diff --git a/runIterativeSolver.m b/runIterativeSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..ce2bcf8b6d8cd180066937d7c5843e07e05907f6
--- /dev/null
+++ b/runIterativeSolver.m
@@ -0,0 +1,145 @@
+function leveldata = runIterativeSolver(maxNiter)
+
+ % Proceed input
+ if nargin < 1
+ maxNiter = 1e2;
+ end
+
+ % Polynomial degree
+ p = 4;
+
+ % Number of refinement levels
+ nLevels = 6;
+
+ % Load mesh
+ mesh = Mesh.loadFromGeometry('unitsquare');
+ mesh.refineUniform();
+
+ % Create FE space
+ fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', ':');
+ u = FeFunction(fes);
+ ustar = FeFunction(fes);
+
+ % Setup variational problem
+ blf = BilinearForm();
+ lf = LinearForm();
+ blf.a = Constant(mesh, 1);
+ blf.qra = QuadratureRule.ofOrder(max(2*(p-1), 1));
+ lf.f = Constant(mesh, 1);
+
+ % Choose iterative solver
+ solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+ % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+ % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
+
+ % Initialize LevelData
+ leveldata = LevelData('results');
+ leveldata.problem = 'Poisson';
+ leveldata.domain = 'UnitSquare';
+ leveldata.method = 'Sp_PCG_AdditiveSchwarz';
+
+ for k = 1:nLevels-1
+ % Output
+ fprintf('Assemble data for level = %d\n', k);
+ % Assemble matrices on intermediate meshes
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
+ freeDofs = getFreeDofs(fes);
+ solver.setupSystemMatrix(A(freeDofs,freeDofs));
+ solver.setupRhs(F(freeDofs)); % initializes x0 as well
+ % Mesh refinement
+ mesh.refineUniform();
+ end
+
+ fprintf('Assemble data for level = %d\n', nLevels);
+
+ % Assemble matrices on fine mesh
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
+ freeDofs = getFreeDofs(fes);
+ solver.setupSystemMatrix(A(freeDofs,freeDofs));
+ solver.setupRhs(F(freeDofs)); % initializes x0 as well
+
+ % Compute exact solution
+ fprintf('Do exact solve\n');
+ xstar = A(freeDofs,freeDofs) \ F(freeDofs);
+
+ % Refinement loop
+ nIter = 1;
+ while true
+ ndof = getDofs(fes).nDofs;
+
+ % Compute energy norm
+ errorU = sqrt((xstar - solver.x)' * A(freeDofs,freeDofs) * (xstar - solver.x));
+
+ % Compute contraction factor for energy norm
+ if nIter == 1
+ contraction = Inf;
+ else
+ contraction = errorU / errorUPrevious;
+ end
+
+ % Store variables
+ leveldata.append('ndof', uint32(ndof), ...
+ 'errorU', errorU);
+ leveldata.setAbsolute(nIter, 'nIter', uint32(nIter), ...
+ 'contraction', contraction);
+
+ % Print level information
+ leveldata.printLevel();
+
+ % Plot solution
+ % if p == 2
+ % x = zeros(mesh.nCoordinates + mesh.nEdges, 1);
+ % x(freeDofs) = solver.x;
+ % figure(1);
+ % plotS2(mesh.clone(), x);
+ %
+ % figure(2);
+ % x(freeDofs) = xstar;
+ % plotS2(mesh.clone(), x);
+ % title("Exact algebraic solution")
+ % else
+ % figure(1);
+ % u.setFreeData(solver.x);
+ % plot(u);
+ %
+ % figure(2);
+ % ustar.setFreeData(xstar);
+ % plot(ustar);
+ % end
+ % pause(2)
+
+ % Break condition
+ % if ndof > maxNdof
+ % break;
+ % end
+
+ % Mesh refinement
+ % mesh.refineUniform();
+
+ if nIter > maxNiter
+ break;
+ end
+
+ % One iterative solver step
+ solver.step();
+
+ % Update iteration counter
+ nIter = nIter + 1;
+ errorUPrevious = errorU;
+ end
+
+ % Plot solution
+ u.setFreeData(solver.x);
+ plot(u);
+
+ % Plot error
+ figure();
+ leveldata.plot('nIter');
+
+ % Plot contraction
+ figure();
+ leveldata.plotAbsolute('nIter');
+end
diff --git a/test.m b/test.m
new file mode 100644
index 0000000000000000000000000000000000000000..505e85f64341a4cb0f5b3508fa0e7e2d6eb4d257
--- /dev/null
+++ b/test.m
@@ -0,0 +1,61 @@
+mesh = Mesh.loadFromGeometry('unitsquare');
+%mesh.refineUniform();
+
+fes1 = FeSpace(mesh, HigherOrderH1Fe(1), 'dirichlet', ':');
+fes2 = FeSpace(mesh, HigherOrderH1Fe(15), 'dirichlet', ':');
+
+testMatrix = ones(getDofs(fes2).nDofs, 1);
+inclusionMatrix = MeshProlongation(fes1, fes2);
+inclusionMatrix = full(inclusionMatrix);
+Vertices = 1:getDofs(fes1).nDofs;
+inclusionMatrix2 = interpolateData(eye(length(Vertices)), fes1, fes2)';
+diff = inclusionMatrix - inclusionMatrix2;
+
+function inclusionMatrix = MeshProlongation(fromFes, toFes)
+ %Create the prolongation matrix on the unit triangle
+ unittriangle = Mesh.loadFromGeometry('unittriangle');
+ fromFesUnitTriangle = FeSpace(unittriangle, ...
+ HigherOrderH1Fe(fromFes.finiteElement.order));
+ toFesUnitTriangle = FeSpace(unittriangle, ...
+ HigherOrderH1Fe(toFes.finiteElement.order));
+ Vertices = 1:getDofs(fromFesUnitTriangle).nDofs;
+ unitTriangleInclusionMatrix = ...
+ permute(interpolateData(eye(length(Vertices)),...
+ fromFesUnitTriangle, toFesUnitTriangle), [2 1]);
+ % Get local to global map in unit triangle
+
+ unitTriangleInclusionMatrix = unitTriangleInclusionMatrix( ...
+ getDofs(fromFesUnitTriangle).element2Dofs, ...
+ getDofs(toFesUnitTriangle).element2Dofs);
+ % Get dofs of the finite element spaces on the meshes
+ fromFesDofs = getDofs(fromFes);
+ toFesDofs = getDofs(toFes);
+
+ % Copy the matrix from the unit triangle for each element
+ mat = repmat(unitTriangleInclusionMatrix, ...
+ [1, 1, fromFes.mesh.nElements]);
+
+ % Create index sets for the accumarray
+ temp = fromFesDofs.element2Dofs(:, :);
+ temp2 = toFesDofs.element2Dofs(:, :);
+ I = repmat(permute(temp, [1 3 2]), [1, size(temp2, 1), 1]);
+ J = repmat(permute(temp2, [3 1 2]), [size(temp, 1), 1, 1]);
+ ind = [I(:), J(:)];
+
+ inclusionMatrix = ...
+ accumarray(ind, mat(:), [], @mean, [], true);
+end
+
+function interpolatedData = interpolateData(data, fromFes, toFes)
+ Dofs = 1:getDofs(toFes).nDofs;
+ nComponents = size(data, 2);
+ interpolatedData = zeros(numel(Dofs), nComponents);
+ feFunctionWrapper = FeFunction(fromFes);
+ for k = 1:nComponents
+ feFunctionWrapper.setData(data(:,k));
+ wholeData = nodalInterpolation(feFunctionWrapper, toFes);
+ interpolatedData(:,k) = wholeData(Dofs);
+ end
+end
+
+
diff --git a/tests/TestCoordinateTransformations.m b/tests/TestCoordinateTransformations.m
index 7ea0d17488000c475c1ecddb1d19743c836b0c8e..ff492ae1a221871385f4611d06174e88bf276c27 100644
--- a/tests/TestCoordinateTransformations.m
+++ b/tests/TestCoordinateTransformations.m
@@ -6,7 +6,7 @@ end
methods (TestMethodSetup)
function createUnitTriangleCoordinates(testCase)
- mesh = Mesh.loadFromGeometry('unittriangle');
+ mesh = Mesh.unitTriangle();
testCase.coordinates = mesh.coordinates;
testCase.fatalAssertNumElements(testCase.coordinates, 6);
end
diff --git a/tests/TestFeSpace.m b/tests/TestFeSpace.m
index 579d6eba5d07c1d287a3501aba60b97299219fd3..fa8f58e1c70e52025fb4f913b77c64600ab4f40a 100644
--- a/tests/TestFeSpace.m
+++ b/tests/TestFeSpace.m
@@ -56,7 +56,7 @@ methods (Test, ParameterCombination='sequential')
end
function dofDataIsComputedCorrectly(testCase, fe, nDof, nFreeDof)
- mesh = Mesh.loadFromGeometry('unittriangle');
+ mesh = Mesh.unitTriangle();
fes = FeSpace(mesh, fe);
testCase.verifyEqual(getDofs(fes).nDofs, nDof)
testCase.verifyEqual(nnz(getFreeDofs(fes)), nFreeDof)
diff --git a/tests/TestFiniteElementFunctions.m b/tests/TestFiniteElementFunctions.m
index f208227e5759ddeef015997f8a24a3eea1b81e43..246c94d328c2fe3320052dc57bbf339f4b92cdcf 100644
--- a/tests/TestFiniteElementFunctions.m
+++ b/tests/TestFiniteElementFunctions.m
@@ -81,7 +81,7 @@ methods (Test)
fes = FeSpace(mesh, testCase.u.fes.finiteElement);
qr = QuadratureRule.ofOrder(max(fes.finiteElement.order,1));
v = FeFunction(fes);
- P = FeProlongation(fes);
+ P = MeshProlongation(fes);
for k = 1:min(5, size(getDofs(fes).element2Dofs, 1))
testCase.setToElementwiseBasisFunction(v, k);
before = sum(integrateElement(v, qr));
@@ -96,12 +96,7 @@ end
methods (Access=private)
function refineAndInterpolate(~, mesh, fes, v, val)
v.setData(val);
- if isa(fes.finiteElement, 'LowestOrderH1Fe') ...
- || isa(fes.finiteElement, 'LowestOrderL2Fe')
- P = LoFeProlongation(fes);
- else
- P = FeProlongation(fes);
- end
+ P = Prolongation.chooseFor(fes);
mesh.refineUniform();
v.setData(prolongate(P, v));
end
diff --git a/tests/TestIterativeSolver.m b/tests/TestIterativeSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..6a32e6af2f72e24d59a80be2ceb5d94690187682
--- /dev/null
+++ b/tests/TestIterativeSolver.m
@@ -0,0 +1,94 @@
+classdef TestIterativeSolver < matlab.unittest.TestCase
+
+properties (TestParameter)
+ p = struct('low', 1, 'medium', 2, 'high', 5);
+ variant = struct('CG', ["cg", ""],...
+ 'direct', ["direct", ""], ...
+ 'jacobiPCG', ["pcg", "jacobi"], ...
+ 'iChol', ["pcg", "iChol"], ...
+ 'loAdditiveSchwarzPCG', ["pcg", "additiveSchwarzLowOrder"], ...
+ 'hoAdditiveSchwarzPCG', ["pcg", "additiveSchwarzHighOrder"], ...
+ 'lowMG', ["multigrid", "lowOrderVcycle"], ...
+ 'highMG', ["multigrid", "highOrderVcycle"]);
+end
+
+methods (Test)
+ function firstStepDecreasesNorm(testCase, p, variant)
+ [~, fes, blf, lf] = setupProblem(testCase, p);
+ s = variant(1); v = variant(2);
+ solver = IterativeSolver.choose(fes, blf, s, v);
+
+ [xstar, A, F] = assembleData(testCase, blf, lf, fes);
+ solver.setupSystemMatrix(A);
+ solver.setupRhs(F);
+
+ normBefore = sqrt((xstar - solver.x)' * A * (xstar - solver.x));
+ solver.step();
+ normAfterwards = sqrt((xstar - solver.x)' * A * (xstar - solver.x));
+
+ testCase.verifyGreaterThan(normBefore, normAfterwards);
+ end
+
+ function laterStepDecreasesNorm(testCase, p, variant)
+ [mesh, fes, blf, lf] = setupProblem(testCase, p);
+ s = variant(1); v = variant(2);
+ solver = IterativeSolver.choose(fes, blf, s, v);
+
+ for k = 1:3
+ [xstar, A, F] = assembleData(testCase, blf, lf, fes);
+ solver.setupSystemMatrix(A);
+ if k < 3, mesh.refineUniform(); end
+ end
+
+ solver.setupRhs(F);
+ normBefore = sqrt((xstar - solver.x)' * A * (xstar - solver.x));
+ solver.step();
+ normAfterwards = sqrt((xstar - solver.x)' * A * (xstar - solver.x));
+
+ testCase.verifyGreaterThan(normBefore, normAfterwards);
+ end
+
+ function solverIsLinear(testCase, p, variant)
+ [mesh, fes, blf, lf] = setupProblem(testCase, p);
+ s = variant(1); v = variant(2);
+ solver = IterativeSolver.choose(fes, blf, s, v);
+
+ for k = 1:3
+ [xstar, A, F] = assembleData(testCase, blf, lf, fes);
+ solver.setupSystemMatrix(A);
+ if k < 3, mesh.refineLocally(k); end
+ end
+
+ solver.setupRhs([F, -pi*F], 0*[F, F]);
+ solver.solve();
+
+ testCase.verifyEqual(-pi*solver.x(:,1), solver.x(:,2), 'RelTol', 2e-5);
+ testCase.verifyEqual(solver.x(:,1), xstar, 'RelTol', 2e-5);
+ end
+end
+
+methods (Access=private)
+ function [mesh, fes, blf, lf] = setupProblem(~, p)
+ mesh = Mesh.loadFromGeometry('unitsquare');
+ fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', ':');
+ blf = BilinearForm();
+ lf = LinearForm();
+ blf.a = Constant(mesh, 1);
+ blf.qra = QuadratureRule.ofOrder(max(2*(p-1), 1));
+ blf.c = Constant(mesh, 1);
+ blf.qrc = QuadratureRule.ofOrder(2*p);
+ lf.f = Constant(mesh, 1);
+ end
+
+ function [xstar, A, F] = assembleData(~, blf, lf, fes)
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
+ freeDofs = getFreeDofs(fes);
+ A = A(freeDofs,freeDofs);
+ F = F(freeDofs);
+ xstar = A \ F;
+ end
+end
+
+
+end
diff --git a/tests/TestPatchwiseMatrix.m b/tests/TestPatchwiseMatrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..d44b7b64dd3a32b36406483f4e93764cc8bdd63b
--- /dev/null
+++ b/tests/TestPatchwiseMatrix.m
@@ -0,0 +1,69 @@
+classdef TestPatchwiseMatrix < matlab.unittest.TestCase
+
+properties (TestParameter)
+ p = {1, 2, 3, 4, 5};
+end
+
+methods (Test)
+ function lowestOrderIsDiagonalInverse(testCase)
+ mesh = Mesh.loadFromGeometry('Lshape');
+ fes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', [], 'neumann', ':');
+ blf = BilinearForm();
+ blf.a = Constant(mesh, 1);
+ blf.c = Constant(mesh, 1);
+ for ell = 1:5
+ mesh.refineUniform();
+ A = assemble(blf, fes);
+ B = assemblePatchwise(blf, fes);
+ actual = B \ full(diag(A));
+ testCase.verifyEqual(actual, ones(getDofs(fes).nDofs, 1), 'RelTol', 10*eps)
+ end
+ end
+
+ function dofsAreInCorrectlyManyPatches(testCase, p)
+ mesh = Mesh.loadFromGeometry('Lshape');
+ for ell = 1:5, mesh.refineLocally(1); end
+ fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', 1, 'neumann', 2);
+ blf = BilinearForm();
+ blf.a = Constant(mesh, 1);
+ blf.qra = QuadratureRule.ofOrder(max(2*(p-1),1));
+ A = assemblePatchwise(blf, fes);
+ countOccurrence = accumarray(getPatchDofs(A, ':'), 1, [getDofs(fes).nDofs, 1]);
+
+ isFree = false(getDofs(fes).nDofs, 1);
+ isFree(getFreeDofs(fes)) = true;
+ % vertex dofs should occur exactly once (if free)
+ vertexDofs = createVertexDofs(fes, ':');
+ testCase.verifyEqual(countOccurrence(vertexDofs), 1*isFree(vertexDofs));
+ % edge dofs should occur exactly twice (if free)
+ edgeDofs = createInnerEdgeDofs(fes, ':');
+ testCase.verifyEqual(countOccurrence(edgeDofs), 2*isFree(edgeDofs));
+ % inner dofs should occur exactly three times
+ elementDofs = createInnerElementDofs(fes, ':');
+ testCase.verifyEqual(countOccurrence(elementDofs), 3*isFree(elementDofs));
+ end
+
+ function patchwiseRelatesToFullMatrix(testCase, p)
+ mesh = Mesh.loadFromGeometry('unitsquare');
+ for ell = 1:5, mesh.refineLocally(1); end
+ fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', [], 'neumann', ':');
+ blf = BilinearForm();
+ blf.a = Constant(mesh, 1);
+ blf.qra = QuadratureRule.ofOrder(max(2*(p-1),1));
+ blf.c = MeshFunction(mesh, @(x) 1+x(1,:,:));
+ blf.qrc = QuadratureRule.ofOrder(2*p+1);
+ A = assemble(blf, fes);
+ B = assemblePatchwise(blf, fes);
+ for k = 1:10:51
+ vertexDof = createVertexDofs(fes, k);
+ idx = getPatchDofs(B, k);
+ expected = zeros(getDofs(fes).nDofs, 1);
+ expected(vertexDof) = 1;
+ actual = (B \ expected);
+ actual(idx) = A(idx,idx) * actual(idx);
+ testCase.verifyEqual(actual, expected, 'AbsTol', sqrt(eps))
+ end
+ end
+end
+
+end
diff --git a/tmp/plotS2.m b/tmp/plotS2.m
new file mode 100644
index 0000000000000000000000000000000000000000..92ee44046520096223c78bb8688872588c419297
--- /dev/null
+++ b/tmp/plotS2.m
@@ -0,0 +1,49 @@
+function ax = plotS2(mesh, x)
+%%PLOTS2 plots a piecewise P2 and globally continuous function (Lagrange
+%finite element function) given by the coefficient vector x on the
+%triangulation meshdata
+% ax = PLOTS2(meshdata, x)
+
+% Copyright 2023 Philipp Bringmann
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+
+
+ % Create axes object
+ ax = gca;
+
+ % Uniform red-refinement of triangulation
+ mesh.refineUniform();
+
+ % S1 plot of nodal interpolation on red-refinement
+ if size(x, 2) == 1
+ % Plot of scalar function
+ pat = trisurf(mesh.elements', mesh.coordinates(1,:)',...
+ mesh.coordinates(2,:)', x);
+ if(mesh.nElements > 2000)
+ set(pat, 'EdgeColor', 'none');
+ end
+ elseif size(x, 2) == 2
+ % Plot of vector field
+ set(ax, 'XLim', 2.1*[min(mesh.coordinates(1,:)), max(mesh.coordinates(1,:))],...
+ 'YLim', 1.1*[min(mesh.coordinates(2,:)), max(mesh.coordinates(2,:))]);
+ quiver(mesh.coordinates(1,:), mesh.coordinates(2,:), x(:,1), x(:,2));
+ else
+ error('Invalid size of coefficient vector');
+ end
+
+ % Draw title
+ title(ax, {'S2 function plot'; [num2str(mesh.nEdges + mesh.nCoordinates), ' dofs']});
+end