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..20caad536454a8d708221bb46d73ebd5af66d027 100644
--- a/README.md
+++ b/README.md
@@ -23,13 +23,14 @@ Advantages:
- Easy to use and modify
- Covers a wide range of equations
- Efficient implementation of general polynomial orders
+- Integrated algebraic multilevel solvers
- A lot of built-in convenience functions
## Installation
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
@@ -55,13 +56,15 @@ runtests('tests/', ['IncludeSubfolders', true,] 'ReportCoverageFor', 'lib/')
## Citation
-If you use the MooAFEM package in your own research, please acknowledge this by citing the accompanying [documentation paper](https://arxiv.org/abs/2203.01845):
+If you use the MooAFEM package in your own research, please acknowledge this by citing the accompanying [documentation paper](https://doi.org/10.1016/j.amc.2022.127731):
```
@Article{MooAFEM,
- author = {Michael Innerberger and Dirk Praetorius},
- title = {{MooAFEM}: An object oriented {Matlab} code for higher-order (nonlinear) adaptive {FEM}},
- year = {2022},
- archiveprefix = {arXiv},
- eprint = {2203.01845},
+ author = {Michael Innerberger and Dirk Praetorius},
+ journal = {Applied Mathematics and Computation},
+ title = {{MooAFEM}: An object oriented {Matlab} code for higher-order adaptive {FEM} for (nonlinear) elliptic {PDEs}},
+ year = {2023},
+ pages = {127731},
+ volume = {442},
+ doi = {10.1016/j.amc.2022.127731},
}
-```
+```
diff --git a/examples/adaptiveFEM.m b/examples/adaptiveFEM.m
index 6b895ca67dc9085932990f27b67cb35d4d35504d..ee88317eff1e9bb54138beafd2267428c0bb0b4d 100644
--- a/examples/adaptiveFEM.m
+++ b/examples/adaptiveFEM.m
@@ -14,8 +14,8 @@ midpoint = [1;1;1]/3;
u = FeFunction(fes);
%% problem data
-blf = BilinearForm(fes);
-lf = LinearForm(fes);
+blf = BilinearForm();
+lf = LinearForm();
blf.a = Constant(mesh, [1;0;0;2]);
blf.b = Constant(mesh, [1;1]);
blf.c = MeshFunction(mesh, @(x) sum(x.^2, Dim.Vector));
@@ -25,11 +25,12 @@ w = FeFunction(ncFes);
lf.fvec = CompositeFunction(@(w) [w;-w], w);
%% adaptive loop
-while 1
+meshSufficientlyFine = false;
+while ~meshSufficientlyFine
%% compute data for rhs coefficient and assemble forms
w.setData(eval(cutoff, midpoint));
- A = assemble(blf);
- F = assemble(lf);
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
%% solve FEM system
freeDofs = getFreeDofs(fes);
@@ -40,13 +41,13 @@ while 1
printLogMessage('number of elements: %d, estimator: %.2e', mesh.nElements, sqrt(sum(eta2)));
%% stoping criterion
- if mesh.nElements > nEmax
- break
- end
+ meshSufficientlyFine = (mesh.nElements > nEmax);
%% refine mesh
- marked = markDoerflerBinning(eta2, theta);
- mesh.refineLocally(marked, 'NVB');
+ if ~meshSufficientlyFine
+ marked = markDoerflerBinning(eta2, theta);
+ mesh.refineLocally(marked, 'NVB');
+ end
end
%% plot solution and estimator
@@ -63,7 +64,8 @@ end
% \eta(T)^2 = h_T^2 * || -div(a*Du) + b*Du + c*u - f + div(fvec) ||_{L^2(T)}^2
% + h_T * || [ a*Du - fvec ] ||_{L^2(E)}^2
function indicators = estimate(blf, lf, u)
- mesh = blf.fes.mesh;
+ fes = u.fes;
+ mesh = fes.mesh;
%% compute volume residual element-wise
% Here, one can optimize the estimator for the given situation. In the case
@@ -83,7 +85,7 @@ function indicators = estimate(blf, lf, u)
blf.a, lf.fvec, Gradient(u));
qrEdge = QuadratureRule.ofOrder(1, '1D');
edgeRes = integrateNormalJump(f, qrEdge, @(j) j.^2, {}, ':');
- dirichlet = unique(vertcat(mesh.boundaries{:}));
+ dirichlet = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
edgeRes(dirichlet) = 0;
%% combine the resdiuals suitably
diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..9e6211ebd724ca7f80ce5d711e49ab5a679322c5
--- /dev/null
+++ b/examples/algebraicSolver.m
@@ -0,0 +1,142 @@
+% ******************************************************************************
+% 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);
+ lf.f = Constant(mesh, 1);
+
+ %% set up solver and operator for nested iteration
+ % choose the iterative solver of choice: multigrid (with the variants
+ % lowOrderVcycle and highOrderVcycle), pcg (with jacobi and additive
+ % Schwarz preconditioner), and the cg solver
+ %[solver, P] = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+ %[solver, P] = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+ [solver, P] = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
+ solver.tol = 1e-6;
+ solver.maxIter = 100;
+
+ %% adaptive loop
+ ell = 0;
+ meshSufficientlyFine = false;
+ while ~meshSufficientlyFine
+ %% setup
+ tic;
+ ell = ell + 1;
+ A = assemble(blf, fes);
+ rhs = assemble(lf, fes);
+
+ freeDofs = getFreeDofs(fes);
+ A = A(freeDofs,freeDofs);
+ solver.setupSystemMatrix(A);
+ solver.setupRhs(rhs(freeDofs), u.data(freeDofs)');
+
+ % exact solution as reference
+ tic;
+ x_ex = A \ rhs(freeDofs);
+ directSolveTime(p, ell) = toc;
+
+ %% iterative solve
+ % solver.solve() implements the algebraic solver loop and
+ % terminates automatically once the wanted accuracy is reached
+ % Another option is to manually implement a loop and implement one
+ % solver step with solver.step() and stop by setting a stopping
+ % criterion
+ solverIsConverged = false;
+ while ~solverIsConverged
+ solver.step();
+ u.setFreeData(solver.x)
+ 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(~, 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/convergenceRates.m b/examples/convergenceRates.m
index 2e200e793e389b7501548882e6a5b6b7bbb99103..d054129fbbc44f028d3f3ff90f5f06e06f5e69e8 100644
--- a/examples/convergenceRates.m
+++ b/examples/convergenceRates.m
@@ -14,49 +14,49 @@ for p = 1:pmax
%% setup geometry & spaces
printLogMessage('*** p = %d (of %d) ***', p, pmax)
mesh = Mesh.loadFromGeometry('Lshape');
- fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', 1);
+ fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', 1, 'neumann', 2);
u = FeFunction(fes);
uex = FeFunction(fes);
%% set problem data for -\Delta u = 1 on L-shape
- blf = BilinearForm(fes);
- lf = LinearForm(fes);
+ blf = BilinearForm();
+ lf = LinearForm();
+ % quadrature rules are set automatically based on order of FeSpace
blf.a = Constant(mesh, 1);
- blf.qra = QuadratureRule.ofOrder(max(2*p-2, 1));
-
+ % ... but can also be set manually if required
lf.neumann = MeshFunction(mesh, @exactSolutionNeumannData);
lf.qrNeumann = QuadratureRule.ofOrder(2*p, '1D');
- lf.bndNeumann = 2;
%% adaptive loop
- i = 1;
- while 1
+ ell = 1;
+ meshSufficientlyFine = false;
+ while ~meshSufficientlyFine
%% assemble & solve FEM system
- A = assemble(blf);
- F = assemble(lf);
+ ell = ell + 1;
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
free = getFreeDofs(fes);
u.setFreeData(A(free,free) \ F(free));
uex.setData(nodalInterpolation(MeshFunction(mesh, @exactSolution), fes));
%% estimate error and store data
eta2 = estimate(blf, lf, u);
- nDofs(p,i) = getDofs(fes).nDofs;
- errEst(p,i) = sqrt(sum(eta2));
+ nDofs(p,ell) = getDofs(fes).nDofs;
+ errEst(p,ell) = sqrt(sum(eta2));
deltaU = u.data - uex.data;
- h1Err(p,i) = sqrt(deltaU * A * deltaU');
- nElem(p,i) = mesh.nElements;
- printLogMessage('number of dofs: %d, estimator: %.2e', nDofs(p,i), errEst(p,i));
+ h1Err(p,ell) = sqrt(deltaU * A * deltaU');
+ nElem(p,ell) = mesh.nElements;
+ printLogMessage('number of dofs: %d, estimator: %.2e', nDofs(p,ell), errEst(p,ell));
%% stoping criterion
- if nDofs(p,i) > nDofsMax
- break
- end
+ meshSufficientlyFine = (nDofs(p,ell) > nDofsMax);
%% refine mesh
- marked = markDoerflerSorting(eta2, theta);
- mesh.refineLocally(marked, 'NVB');
- i = i+1;
+ if ~meshSufficientlyFine
+ marked = markDoerflerSorting(eta2, theta);
+ mesh.refineLocally(marked, 'NVB');
+ end
end
end
@@ -87,9 +87,10 @@ end
% \eta(T)^2 = h_T^2 * || \Delta u + f ||_{L^2(T)}^2
% + h_T * || [[Du * n]] ||_{L^2(E) \cap \Omega}^2
% + h_T * || Du * n - phi ||_{L^2(E) \cap \Gamma_N}^2
-function indicators = estimate(blf, lf, u)
- p = blf.fes.finiteElement.order;
- mesh = blf.fes.mesh;
+function indicators = estimate(~, lf, u)
+ fes = u.fes;
+ p = fes.finiteElement.order;
+ mesh = fes.mesh;
% compute volume residual element-wise
% For p=1, the diffusion term vanishes in the residual.
@@ -103,9 +104,11 @@ function indicators = estimate(blf, lf, u)
% compute edge residual edge-wise
qr = QuadratureRule.ofOrder(p, '1D');
+ dirichletBnd = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
+ neumannBnd = getCombinedBndEdges(mesh, fes.bnd.neumann);
edgeRes = integrateNormalJump(Gradient(u), qr, ...
- @(j) zeros(size(j)), {}, mesh.boundaries{1}, ... % no jump on dirichlet edges
- @(j,phi) j-phi, {lf.neumann}, mesh.boundaries{2},...% neumann jump on neumann edges
+ @(j) zeros(size(j)), {}, dirichletBnd, ... % no jump on dirichlet edges
+ @(j,phi) j-phi, {lf.neumann}, neumannBnd,... % neumann jump on neumann edges
@(j) j.^2, {}, ':'); % normal jump on inner edges
% combine the resdiuals suitably
diff --git a/examples/femSystemAssembly.m b/examples/femSystemAssembly.m
index 34490af04286e65ba53bc6491b4e388ca132956c..7650b1d158effae48f344d98dd38de79034234e5 100644
--- a/examples/femSystemAssembly.m
+++ b/examples/femSystemAssembly.m
@@ -11,8 +11,8 @@ fes = FeSpace(mesh, LowestOrderH1Fe());
% One can define a bilinear form 'A' and a linear form 'F', such that
% A(u,v) = \int a*Du*Dv + b*Du*v + c*u*v dx,
% F(v) = \int f*v + fvec*Dv dx.
-A = BilinearForm(fes);
-F = LinearForm(fes);
+A = BilinearForm();
+F = LinearForm();
%% coefficients
% The coefficients a, b, c and f, fvec can be arbitrary Evaluables.
@@ -22,8 +22,8 @@ F.f = Constant(mesh, 1);
%% assembly and solution
% Once the coefficients are fixed, the bilinear form and linear form can be
% assembled to yield a (sparse) matrix and a vector, respectively.
-Amat = assemble(A);
-Fvec = assemble(F);
+Amat = assemble(A, fes);
+Fvec = assemble(F, fes);
% This allows for easy solution of the corresponding linear system.
freeDofs = getFreeDofs(fes);
@@ -39,8 +39,8 @@ A.c = CompositeFunction(@(x) x.^2, u);
F.f = [];
F.fvec = MeshFunction(mesh, @(x) -(0.5-x).^2);
-Amat = assemble(A);
-Fvec = assemble(F);
+Amat = assemble(A, fes);
+Fvec = assemble(F, fes);
w = FeFunction(fes);
w.setFreeData( Amat(freeDofs, freeDofs) \ Fvec(freeDofs) );
@@ -49,8 +49,8 @@ plot(w)
%% other FeSpaces
% Of course, also other FeSpaces can be used to assemble the forms.
fes2 = FeSpace(mesh, LowestOrderL2Fe());
-M = BilinearForm(fes2);
-G = LinearForm(fes2);
+M = BilinearForm();
+G = LinearForm();
% This realizes mass matrix and integration for piecewise constant functions.
% NOTE: All assembly operations are implemented in terms of local shape
@@ -59,8 +59,8 @@ G = LinearForm(fes2);
M.c = Constant(mesh, 1);
G.f = Constant(mesh, 1);
-Mmat = assemble(M);
-Gvec = assemble(G);
+Mmat = assemble(M, fes2);
+Gvec = assemble(G, fes2);
v = FeFunction(fes2);
v.setData(ones(mesh.nElements,1));
diff --git a/examples/goafemIterativeSolver.m b/examples/goafemIterativeSolver.m
new file mode 100644
index 0000000000000000000000000000000000000000..aed69d04e8ac2983c895c06b27299ea46a611128
--- /dev/null
+++ b/examples/goafemIterativeSolver.m
@@ -0,0 +1,149 @@
+% ******************************************************************************
+% 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);
+
+ 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(~, 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 e966ffa0bc7ff92cf2f609649d2436c66ea0e58c..a07806216b2b6a115b00a751585e3797c48ba33d 100644
--- a/examples/goalOrientedAFEM.m
+++ b/examples/goalOrientedAFEM.m
@@ -22,31 +22,31 @@ for p = 1:pmax
z = FeFunction(fes);
%% set problem data given in [Mommer, Stevenson; 2009]
- blf = BilinearForm(fes);
+ % quadrature orders are set automatically based on order of FeSpace
+ blf = BilinearForm();
blf.a = Constant(mesh, 1);
- blf.qra = QuadratureRule.ofOrder(max(2*p-2, 1));
chiT1 = MeshFunction(mesh, @(x) sum(x, Dim.Vector) < 1/2);
v.setData(nodalInterpolation(chiT1, ncFes));
- lfF = LinearForm(fes);
+ lfF = LinearForm();
lfF.fvec = CompositeFunction(@(v) [v;zeros(size(v))], v);
- lfF.qrfvec = QuadratureRule.ofOrder(max(p-1, 1));
chiT2 = MeshFunction(mesh, @(x) sum(x, Dim.Vector) > 3/2);
w.setData(nodalInterpolation(chiT2, ncFes));
- lfG = LinearForm(fes);
+ lfG = LinearForm();
lfG.fvec = CompositeFunction(@(w) [-w;zeros(size(w))], w);
- lfG.qrfvec = QuadratureRule.ofOrder(max(p-1, 1));
%% set up lifting operators for rhs FEM-data
- P = LoFeProlongation(ncFes);
+ P = LoMeshProlongation(ncFes);
%% adaptive loop
- i = 1;
- while 1
+ ell = 1;
+ meshSufficientlyFine = false;
+ while ~meshSufficientlyFine
%% assemble & solve FEM system
- A = assemble(blf);
- rhs = [assemble(lfF), assemble(lfG)];
+ ell = ell + 1;
+ A = assemble(blf, fes);
+ rhs = [assemble(lfF, fes), assemble(lfG, fes)];
freeDofs = getFreeDofs(fes);
uz = A(freeDofs,freeDofs) \ rhs(freeDofs,:);
u.setFreeData(uz(:,1));
@@ -55,22 +55,21 @@ for p = 1:pmax
%% estimate error and store data
eta2 = estimate(blf, lfF, u);
zeta2 = estimate(blf, lfG, z);
- nDofs(p,i) = getDofs(fes).nDofs;
- nElem(p,i) = mesh.nElements;
- goalErrEst(p,i) = sqrt(sum(eta2)*sum(zeta2));
- printLogMessage('number of dofs: %d, estimator: %.2e', nDofs(p,i), goalErrEst(p,i));
+ nDofs(p,ell) = getDofs(fes).nDofs;
+ nElem(p,ell) = mesh.nElements;
+ goalErrEst(p,ell) = sqrt(sum(eta2)*sum(zeta2));
+ printLogMessage('number of dofs: %d, estimator: %.2e', nDofs(p,ell), goalErrEst(p,ell));
%% stoping criterion
- if nDofs(p,i) > nDofsMax
- break
- end
+ meshSufficientlyFine = (nDofs(p,ell) > nDofsMax);
%% refine mesh
- marked = markGoafemMS(eta2, zeta2, theta);
- mesh.refineLocally(marked, 'NVB');
- i = i+1;
- v.setData(prolongate(P, v));
- w.setData(prolongate(P, w));
+ if ~meshSufficientlyFine
+ marked = markGoafemMS(eta2, zeta2, theta);
+ mesh.refineLocally(marked, 'NVB');
+ v.setData(prolongate(P, v));
+ w.setData(prolongate(P, w));
+ end
end
end
@@ -94,9 +93,10 @@ title(['goal error estimator over number of dofs'])
%% local function for residual a posteriori error estimation
% \eta(T)^2 = h_T^2 * || \Delta u ||_{L^2(T)}^2
% + h_T * || [[(Du - fvec) * n]] ||_{L^2(E) \cap \Omega}^2
-function indicators = estimate(blf, lf, u)
- p = blf.fes.finiteElement.order;
- mesh = blf.fes.mesh;
+function indicators = estimate(~, lf, u)
+ fes = u.fes;
+ p = fes.finiteElement.order;
+ mesh = fes.mesh;
trafo = getAffineTransformation(mesh);
% compute volume residual element-wise
@@ -112,9 +112,10 @@ function indicators = estimate(blf, lf, u)
% compute edge residual edge-wise
qr = QuadratureRule.ofOrder(max(p-1, 1), '1D');
f = CompositeFunction(@(p,fvec) p-fvec, Gradient(u), lf.fvec);
+ dirichletBnd = getCombinedBndEdges(mesh, fes.bnd.dirichlet);
edgeResidual = integrateNormalJump(f, qr, ...
- @(j) zeros(size(j)), {}, mesh.boundaries{1}, ... % no jump on dirichlet edges
- @(j) j.^2, {}, ':'); % normal jump on inner edges
+ @(j) zeros(size(j)), {}, dirichletBnd, ... % no jump on dirichlet edges
+ @(j) j.^2, {}, ':'); % normal jump on inner edges
% combine the resdiuals suitably
hT = sqrt(trafo.area);
diff --git a/examples/iterativeLinearization.m b/examples/iterativeLinearization.m
index 9c9ea5efaf0902140faec7b93cc08eefcb30ef3e..57a6f1782c0382497796e65eda21b310ba3677e7 100644
--- a/examples/iterativeLinearization.m
+++ b/examples/iterativeLinearization.m
@@ -5,7 +5,7 @@
%% paramters
nElemMax = 1e4;
theta = 0.5;
-lambda = 0.5;
+lambda = 0.1;
deltaZ = 0.1;
linearizations = ["zarantonello", "kacanov", "newton"];
@@ -16,65 +16,64 @@ for k = 1:length(linearizations)
printLogMessage('*** starting %s iteration', linearizations(k))
mesh = Mesh.loadFromGeometry('Lshape');
mesh.refineUniform(2);
- fes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', 'all');
+ 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, linearizations(k));
- eDensity = CompositeFunction(@(p) muIntegral(vectorProduct(p, p)), Du);
+ [blf, lf] = setProblemData(mesh, Du, g, linearizations(k));
+ eDensity = CompositeFunction(@(p) 1/2*muIntegral(vectorProduct(p, p)), Du);
u.setData(0);
%% nested iteration
- P = LoFeProlongation(fes);
+ P = LoMeshProlongation(fes);
%% adaptive algorithm
- i = 1;
+ ell = 0;
Eold = 0;
qr = QuadratureRule.ofOrder(1);
- while 1
+ meshSufficientlyFine = false;
+ while ~meshSufficientlyFine
% 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
+ nIterations = 0;
+ solverIsConverged = false;
+ while ~solverIsConverged
%% solve linearized problem
- A = assemble(blf);
- F = assemble(lf);
+ ell = ell + 1;
+ nIterations = nIterations + 1;
+
+ A = assemble(blf, fes);
+ F = assemble(lf, fes);
updateDataU(u, deltaZ, A, F, linearizations(k));
%% compute remainder of error estimator and energy
edge = estimateEdge(u, mesh);
eta2 = combineEstimators(vol, edge, mesh);
- eta(k,i) = sqrt(sum(eta2));
- nElem(k,i) = mesh.nElements;
-
- if i > 1, Eold = E(k,i-1); end
- E(k,i) = sum(integrateElement(eDensity, qr)) + u.data * F;
+ eta(k,ell) = sqrt(sum(eta2));
+ nElem(k,ell) = mesh.nElements;
+ E(k,ell) = sum(integrateElement(eDensity, qr)) + u.data * F;
%% stopping criterion of the iterative solver
- if sqrt(E(k,i) - Eold) <= lambda * eta(k,i)
- break
- end
- nIterations = nIterations + 1;
- i = i+1;
+ solverIsConverged = (sqrt(E(k,ell) - Eold) <= lambda * eta(k,ell));
+ Eold = E(k,ell);
end
%% stoping criterion for the AFEM algorithm
printLogMessage('number of elements: %d, iterations: %d, estimator: %.2e', ...
- nElem(k,i), nIterations, eta(k,i));
- if nElem(k,i) > nElemMax
- break
- end
+ nElem(k,ell), nIterations, eta(k,ell));
+ meshSufficientlyFine = (nElem(k,ell) > nElemMax);
%% refine mesh
- marked = markDoerflerSorting(eta2, theta);
- mesh.refineLocally(marked, 'NVB');
- u.setData(prolongate(P, u));
- i = i+1;
+ if ~meshSufficientlyFine
+ marked = markDoerflerSorting(eta2, theta);
+ mesh.refineLocally(marked, 'NVB');
+ u.setData(prolongate(P, u));
+ end
end
end
@@ -98,9 +97,9 @@ function plotData(xdata, xlab, ydata, ylab, names)
end
%% local function for problem data
-function [blf, lf] = setProblemData(mesh, fes, Du, g, linearization)
- blf = BilinearForm(fes);
- lf = LinearForm(fes);
+function [blf, lf] = setProblemData(mesh, Du, g, linearization)
+ blf = BilinearForm();
+ lf = LinearForm();
switch linearization
case "zarantonello"
@@ -117,8 +116,6 @@ function [blf, lf] = setProblemData(mesh, fes, Du, g, linearization)
lf.f = g;
lf.fvec = CompositeFunction(@(p) -mu(vectorProduct(p, p)) .* p, Du);
end
-
- [blf.qra, lf.qrf, lf.qrfvec] = deal(QuadratureRule.ofOrder(1));
end
function updateDataU(u, deltaZ, A, F, linearization)
@@ -146,7 +143,7 @@ end
function edge = estimateEdge(u, mesh)
qr = QuadratureRule.ofOrder(1, '1D');
- dirichlet = vertcat(mesh.boundaries{:});
+ dirichlet = getCombinedBndEdges(mesh, u.fes.bnd.dirichlet);
f = CompositeFunction(@(p) mu(vectorProduct(p,p)).*p, Gradient(u));
edge = integrateNormalJump(f, qr, ...
diff --git a/examples/meshGenerationAndPlotting.m b/examples/meshGenerationAndPlotting.m
index cb8595cfa9ff9e512e5b3760ef02e72bb019bfae..97341b0a3e3c8b9185098464113faa8da1c5f6f2 100644
--- a/examples/meshGenerationAndPlotting.m
+++ b/examples/meshGenerationAndPlotting.m
@@ -8,10 +8,8 @@
% constructed easily, given the name of the subfolder.
mesh = Mesh.loadFromGeometry('Lshape');
-% A mesh can also be generated from the data by hand (not recommended, except
-% for copying).
-boundaries = cellfun(@(bnd) mesh.edges(:,bnd), mesh.boundaries, 'UniformOutput', false);
-referenceMesh = Mesh(mesh.coordinates, mesh.elements, boundaries);
+% A mesh can also be cloned.
+referenceMesh = clone(mesh);
%% modification
% The data of a mesh can be transformed (translate, scale, rotate, skew) easily.
diff --git a/examples/nonlinearFEM.m b/examples/nonlinearFEM.m
index 7a93797b4d6c6f2bb4c315960926adf9e788c9bb..1fbbd3f462c45f11ce2f5cebc4fde68415d4599a 100644
--- a/examples/nonlinearFEM.m
+++ b/examples/nonlinearFEM.m
@@ -24,10 +24,10 @@ uexact = MeshFunction(mesh, @(x) prod(sin([pi;2*pi].*x), Dim.Vector));
% which change and those which don't.
% The following iteration is based on the approximation
% (u^{n+1})^3 ~ 3*(u^{n})^2 * u^{n+1} - 2*(u^{n})^3.
-blfFixed = BilinearForm(fes);
-blfVariable = BilinearForm(fes);
-lfFixed = LinearForm(fes);
-lfVariable = LinearForm(fes);
+blfFixed = BilinearForm();
+blfVariable = BilinearForm();
+lfFixed = LinearForm();
+lfVariable = LinearForm();
blfFixed.a = Constant(mesh, 1e-1);
lfFixed.f = CompositeFunction(@(uex) 1e-1*5*pi^2*uex + uex.^3, uexact);
@@ -38,8 +38,8 @@ blfVariable.qrc = QuadratureRule.ofOrder(4);
lfVariable.f = CompositeFunction(@(u) 2*u.^3, u);
lfVariable.qrf = QuadratureRule.ofOrder(4);
-B = assemble(blfFixed);
-G = assemble(lfFixed);
+B = assemble(blfFixed, fes);
+G = assemble(lfFixed, fes);
%% iterative solution of non-linear equation
n = 1;
@@ -47,8 +47,8 @@ while n <= nMax
% The variable part of the (bi-)linear form needs to be assembled in every
% iteration, since the underlying data change. Since FeFunctions are
% handles, the current state of u is taken into account in the assembly.
- A = B + assemble(blfVariable);
- F = G + assemble(lfVariable);
+ A = B + assemble(blfVariable, fes);
+ F = G + assemble(lfVariable, fes);
% If we update the data of u, also coefficients which depend on u are updated.
freeDofs = getFreeDofs(fes);
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/lib/assembly/@BilinearForm/BilinearForm.m b/lib/assembly/@BilinearForm/BilinearForm.m
index fbd4f0e807b96750ef9c6f8bf92fb44a780a6f50..5387cef70ff48a924ebfca72f8ce7e23814520a9 100644
--- a/lib/assembly/@BilinearForm/BilinearForm.m
+++ b/lib/assembly/@BilinearForm/BilinearForm.m
@@ -4,37 +4,25 @@
classdef BilinearForm < handle
%% properties
- properties (GetAccess=public, SetAccess=protected)
- fes
- end
-
properties (Access=public)
a {mustBeEvaluableOrEmpty} = [] % Diffusion matrix (or scalar)
b {mustBeEvaluableOrEmpty} = [] % Convection vector field
c {mustBeEvaluableOrEmpty} = [] % Reaction coefficient
robin {mustBeEvaluableOrEmpty} = [] % Robin coefficient
- qrNeumann (1,1) QuadratureRule = QuadratureRule.ofOrder(1, '1D')
- qra (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
- qrb (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
- qrc (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
- qrRobin (1,1) QuadratureRule = QuadratureRule.ofOrder(1, '1D')
- bndRobin {mustBeIndexVector} = []
+ qra (1,1) QuadratureRule = UnspecifiedQR
+ qrb (1,1) QuadratureRule = UnspecifiedQR
+ qrc (1,1) QuadratureRule = UnspecifiedQR
+ qrRobin (1,1) QuadratureRule = UnspecifiedQR
end
%% methods
- methods
- function obj = BilinearForm(fes)
- % BilinearForm Construct bilinear form from FeSpace.
- arguments
- fes (1,1) FeSpace
- end
-
- obj.fes = fes;
- end
+ methods (Access=public)
+ mat = assemble(obj, fes);
+ mat = assemblePatchwise(obj, fes);
end
- methods (Access=public)
- mat = assemble(obj);
+ 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 6fb18b9d6539655e40110c254103ca20c5b0ddf0..cf56fa2b1c8abeefdadafcc345d9abe79c3b006e 100644
--- a/lib/assembly/@BilinearForm/assemble.m
+++ b/lib/assembly/@BilinearForm/assemble.m
@@ -1,41 +1,17 @@
% assemble Assembles matrix of bilinear form.
%
-% A = assemble(obj) assembles the data of the given BilinearForm and
-% returns a sparse matrix A.
+% A = assemble(obj, fes) assembles the data of the given BilinearForm for a
+% FeSpace fes and returns a sparse matrix A.
%
% See also: BilinearForm
-function mat = assemble(obj)
-%% setup
-if all([isempty(obj.a), isempty(obj.b), isempty(obj.c), isempty(obj.robin)])
- error('BilinearForm:NoCoefficients', 'No coefficients are given.')
+function mat = assemble(obj, fes)
+arguments
+ obj
+ fes FeSpace
end
-fes = obj.fes;
-
-%% 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)
@@ -49,11 +25,13 @@ end
%% integrate and add boundary data
% Robin data
if ~isempty(obj.robin)
+ qrRobin = QuadratureRule.ifEmptyOfOrder(obj.qrRobin, 2*fes.finiteElement.order, "1D");
+ phi = TestFunction(fes);
f = CompositeFunction(@BilinearForm.robinPart, obj.robin, phi);
- idx = getCombinedBndEdges(fes.mesh, obj.bndRobin);
+ idx = getCombinedBndEdges(fes.mesh, fes.bnd.robin);
[I, J] = getLocalDofs(size(dofs.edge2Dofs, Dim.Vector));
mat = mat + sparse(dofs.edge2Dofs(I,idx), dofs.edge2Dofs(J,idx), ...
- integrateEdge(f, obj.qrRobin, idx), dofs.nDofs, dofs.nDofs);
+ integrateEdge(f, qrRobin, idx), dofs.nDofs, dofs.nDofs);
end
end
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..3ff54fe51bae08ee2c498d8b9fbf05dce43c2200
--- /dev/null
+++ b/lib/assembly/@BilinearForm/computeVolumeData.m
@@ -0,0 +1,45 @@
+% 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);
+p = fes.finiteElement.order;
+
+% diffusion
+if ~isempty(obj.a)
+ qr = QuadratureRule.ifEmptyOfOrder(obj.qra, 2*(p-1));
+ f = CompositeFunction(@BilinearForm.diffusionPart, obj.a, Dphi);
+ data = data + integrateElement(f, qr);
+end
+
+% convection
+if ~isempty(obj.b)
+ qr = QuadratureRule.ifEmptyOfOrder(obj.qrb, 2*p-1);
+ f = CompositeFunction(@BilinearForm.convectionPart, obj.b, Dphi, phi);
+ data = data + integrateElement(f, qr);
+end
+
+% reaction
+if ~isempty(obj.c)
+ qr = QuadratureRule.ifEmptyOfOrder(obj.qrc, 2*p);
+ f = CompositeFunction(@BilinearForm.reactionPart, obj.c, phi);
+ data = data + integrateElement(f, qr);
+end
+
+end
diff --git a/lib/assembly/@LinearForm/LinearForm.m b/lib/assembly/@LinearForm/LinearForm.m
index 7fc9ddc830b59b8297d949a37b77d16d9a0d2b9b..02e33ad5de8017998483fb554848153c9d7c9d89 100644
--- a/lib/assembly/@LinearForm/LinearForm.m
+++ b/lib/assembly/@LinearForm/LinearForm.m
@@ -4,37 +4,20 @@
classdef LinearForm < handle
%% properties
- properties (GetAccess=public, SetAccess=protected)
- fes
- end
-
properties (Access=public)
f {mustBeEvaluableOrEmpty} = [] % Volume load
fvec {mustBeEvaluableOrEmpty} = [] % Generalized volume load
robin {mustBeEvaluableOrEmpty} = [] % Robin boundary load
neumann {mustBeEvaluableOrEmpty} = [] % Neumann boundary load
- qrf (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
- qrfvec (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
- qrRobin (1,1) QuadratureRule = QuadratureRule.ofOrder(1, '1D')
- qrNeumann (1,1) QuadratureRule = QuadratureRule.ofOrder(1, '1D')
- bndRobin {mustBeIndexVector} = []
- bndNeumann {mustBeIndexVector} = []
+ qrf (1,1) QuadratureRule = UnspecifiedQR
+ qrfvec (1,1) QuadratureRule = UnspecifiedQR
+ qrRobin (1,1) QuadratureRule = UnspecifiedQR
+ qrNeumann (1,1) QuadratureRule = UnspecifiedQR
end
%% methods
- methods
- function obj = LinearForm(fes)
- % LinearForm Construct bilinear form from FeSpace.
- arguments
- fes (1,1) FeSpace
- end
-
- obj.fes = fes;
- end
- end
-
methods (Access=public)
- vec = assemble(obj);
+ vec = assemble(obj, fes);
end
methods (Static, Access=protected)
diff --git a/lib/assembly/@LinearForm/assemble.m b/lib/assembly/@LinearForm/assemble.m
index 933cd7e7d8e1ee43eba3be8df4429f3b8871dba4..af9f28a0ce2b3ec69a36453318cd67dff2aef1d2 100644
--- a/lib/assembly/@LinearForm/assemble.m
+++ b/lib/assembly/@LinearForm/assemble.m
@@ -1,32 +1,38 @@
% assemble Assembles vector of linear form.
%
-% F = assemble(obj) assembles the data of the given LinearForm and
-% returns a vector F.
+% F = assemble(obj, fes) assembles the data of the given LinearForm for a
+% FeSpace fes and returns a vector F.
%
% See also: LinearForm
-function vec = assemble(obj)
+function vec = assemble(obj, fes)
+arguments
+ obj
+ fes FeSpace
+end
%% setup
if all([isempty(obj.f), isempty(obj.fvec), isempty(obj.neumann), isempty(obj.robin)])
error('LinearForm:NoCoefficients', 'No coefficients are given.')
end
-fes = obj.fes;
%% integrate element data
vec = 0;
phi = TestFunction(fes);
Dphi = TestFunctionGradient(fes);
+p = fes.finiteElement.order;
% scalar rhs
if ~isempty(obj.f)
+ qr = QuadratureRule.ifEmptyOfOrder(obj.qrf, p);
f = CompositeFunction(@LinearForm.scalarPart, obj.f, phi);
- vec = vec + integrateElement(f, obj.qrf);
+ vec = vec + integrateElement(f, qr);
end
% vector rhs
if ~isempty(obj.fvec)
+ qr = QuadratureRule.ifEmptyOfOrder(obj.qrfvec, p-1);
f = CompositeFunction(@LinearForm.vectorPart, obj.fvec, Dphi);
- vec = vec + integrateElement(f, obj.qrfvec);
+ vec = vec + integrateElement(f, qr);
end
%% construct rhs vector from element data
@@ -36,18 +42,20 @@ vec = accumarray(dofs.element2Dofs(:), vec(:), [dofs.nDofs, 1]);
%% integrate and add boundary data
% Robin data
if ~isempty(obj.robin)
+ qr = QuadratureRule.ifEmptyOfOrder(obj.qrRobin, p, "1D");
f = CompositeFunction(@LinearForm.robinPart, obj.robin, phi);
- idx = getCombinedBndEdges(fes.mesh, obj.bndRobin);
- edgeData = integrateEdge(f, obj.qrRobin, idx);
+ idx = getCombinedBndEdges(fes.mesh, fes.bnd.robin);
+ edgeData = integrateEdge(f, qr, idx);
edge2Dofs = dofs.edge2Dofs(:,idx);
vec = vec + accumarray(edge2Dofs(:), edgeData(:), [dofs.nDofs, 1]);
end
% Neumann data
if ~isempty(obj.neumann)
+ qr = QuadratureRule.ifEmptyOfOrder(obj.qrNeumann, p-1, "1D");
f = CompositeFunction(@LinearForm.neumannPart, obj.neumann, phi);
- idx = getCombinedBndEdges(fes.mesh, obj.bndNeumann);
- edgeData = integrateEdge(f, obj.qrNeumann, idx);
+ idx = getCombinedBndEdges(fes.mesh, fes.bnd.neumann);
+ edgeData = integrateEdge(f, qr, idx);
edge2Dofs = dofs.edge2Dofs(:,idx);
vec = vec + accumarray(edge2Dofs(:), edgeData(:), [dofs.nDofs, 1]);
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 02333f174279cfbad7b491361d1ec0ab8b35d636..45fe85963d81fdba976e9c5ccac0a76d07898d22 100644
--- a/lib/mesh/@Mesh/Mesh.m
+++ b/lib/mesh/@Mesh/Mesh.m
@@ -33,6 +33,7 @@ classdef Mesh < handle
IsAboutToRefine
JustRefined
RefineCompleted
+ CompletelyChanged
end
%% public methods
@@ -68,6 +69,8 @@ classdef Mesh < handle
refineUniform(obj, n, method)
changeRefinementEdge(obj, newRefinementEdge)
edges = getCombinedBndEdges(obj, idx)
+ clonedMesh = clone(obj);
+ copyData(obj, otherMesh);
end
%% protected methods
@@ -97,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/clone.m b/lib/mesh/@Mesh/clone.m
new file mode 100644
index 0000000000000000000000000000000000000000..8ccdbd4bed335206eaf3d0bface48673b467be23
--- /dev/null
+++ b/lib/mesh/@Mesh/clone.m
@@ -0,0 +1,12 @@
+% clone Create new mesh with same data
+%
+% newMesh = clone(obj)
+%
+% See also Mesh
+
+function newMesh = clone(obj)
+
+newMesh = Mesh([], [], []);
+copyData(newMesh, obj);
+
+end
diff --git a/lib/mesh/@Mesh/copyData.m b/lib/mesh/@Mesh/copyData.m
new file mode 100644
index 0000000000000000000000000000000000000000..f3dc190add9c97eee3bd7fa9561cdf3543538962
--- /dev/null
+++ b/lib/mesh/@Mesh/copyData.m
@@ -0,0 +1,22 @@
+% copyData Copy data from other Mesh object
+%
+% copyData(obj, otherMesh) updates the data stored in the calling mesh based
+% on other mesh.
+%
+% See also Mesh
+
+function copyData(obj, otherMesh)
+
+assert(isa(otherMesh, 'Mesh'))
+
+% if it is a mesh, data must be valid -> just copy
+obj.coordinates = otherMesh.coordinates;
+obj.elements = otherMesh.elements;
+obj.edges = otherMesh.edges;
+obj.element2edges = otherMesh.element2edges;
+obj.flipEdges = otherMesh.flipEdges;
+obj.boundaries = otherMesh.boundaries;
+
+obj.notify('CompletelyChanged')
+
+end
diff --git a/lib/mesh/@Mesh/saveTikzConforming.m b/lib/mesh/@Mesh/saveTikzConforming.m
index dd42bdd3b4e85209e0d3abf1af7c8f2a4668fbd0..7758989a2348e1d4a8e8e31bb4140fa31c5666c7 100644
--- a/lib/mesh/@Mesh/saveTikzConforming.m
+++ b/lib/mesh/@Mesh/saveTikzConforming.m
@@ -21,7 +21,7 @@ end
%% open file
fileID = fopen(fileName, 'a');
assert(fileID ~= -1, 'Cannot open file %s.', fileName)
-fprintf(fileID, '\\tikzset{meshStyle/.style={very thin, black, fill=none}}\n\n');
+fprintf(fileID, '\\tikzset{meshStyle/.style={very thin, black, fill=none, line join=bevel}}\n\n');
%% define coordinates in tikz
fprintf(fileID, '\\coordinate (P%d) at (%.4g,%.4g);\n', [1:obj.nCoordinates; obj.coordinates]);
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/quadrature/quadratureRules/QuadratureRule.m b/lib/quadrature/quadratureRules/QuadratureRule.m
index 38d2340340b1e091e0ce185a3ff0d45e65f630ad..ee69a439465fceca3d349bbd8baa3e8bf3bfe0ec 100644
--- a/lib/quadrature/quadratureRules/QuadratureRule.m
+++ b/lib/quadrature/quadratureRules/QuadratureRule.m
@@ -39,6 +39,7 @@ classdef QuadratureRule
end
methods(Static)
+ % Return quadrature rule of specified order and dimension
function obj = ofOrder(order, dim)
arguments
order {mustBePositive(order)}
@@ -56,6 +57,23 @@ classdef QuadratureRule
obj = QuadratureRule(bary, weights);
end
+
+ % Return given quadrature rule or create new one if given rule is
+ % an UnspecifiedQR
+ function qr = ifEmptyOfOrder(qrOrUnspecified, order, dim)
+ arguments
+ qrOrUnspecified (1,1) QuadratureRule
+ order (1,1) double
+ dim {mustBeTextScalar, mustBeMember(dim, {'1D', '2D'})} = '2D'
+ end
+
+ if isa(qrOrUnspecified, 'UnspecifiedQR')
+ order = max(order, 1);
+ qr = QuadratureRule.ofOrder(order, dim);
+ else
+ qr = qrOrUnspecified;
+ end
+ end
end
end
diff --git a/lib/quadrature/quadratureRules/UnspecifiedQR.m b/lib/quadrature/quadratureRules/UnspecifiedQR.m
new file mode 100644
index 0000000000000000000000000000000000000000..9979bc84b7d7eb2a17ae0a6d0f1591b267c94aea
--- /dev/null
+++ b/lib/quadrature/quadratureRules/UnspecifiedQR.m
@@ -0,0 +1,14 @@
+% UnspecifiedQR Null object as placeholder for not-set quadrature rule.
+
+classdef UnspecifiedQR < QuadratureRule
+ %% methods
+ methods (Access=public)
+ function obj = UnspecifiedQR()
+ obj = obj@QuadratureRule(Barycentric1D([0; 1]), 1);
+ obj.bary = [];
+ obj.weights = [];
+ obj.nNodes = 0;
+ obj.dim = NaN;
+ end
+ end
+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 d0fbc2da0753f1f10c823febdc99020bd0488930..39413fd4187254266dcd80e86aa98cad588a2692 100644
--- a/lib/spaces/@FeSpace/FeSpace.m
+++ b/lib/spaces/@FeSpace/FeSpace.m
@@ -7,7 +7,7 @@ classdef FeSpace < handle
properties (GetAccess=public, SetAccess=protected)
mesh
finiteElement
- dirichlet
+ bnd
end
properties (Access=private)
@@ -22,20 +22,22 @@ classdef FeSpace < handle
%% methods
methods (Access=public)
- function obj = FeSpace(mesh, fe, opt)
+ function obj = FeSpace(mesh, fe, bnd)
% Construct FeSpace from given mesh and finite element.
%
% fes = FeSpace(mesh, fe) constructs FeSpace from FiniteElement fe
% on mesh, whole boundary is dirichlet boundary.
%
- % fes = FeSpace(mesh, fe, 'dirichlet', dir) Dirichlet boundary is
- % given by dir, which is either an array of boundary indices,
- % or 'all'.
+ % fes = FeSpace(mesh, fe, bndName, bndIdx, ...) Boundaries can be
+ % given by tuples consisting of the name ('dirichlet',
+ % 'neumann', 'robin') and an index vector.
arguments
mesh (1,1) Mesh
fe (1,1) FiniteElement
- opt.dirichlet = 'all'
+ bnd.dirichlet {mustBeIndexVector} = []
+ bnd.neumann {mustBeIndexVector} = []
+ bnd.robin {mustBeIndexVector} = []
end
% register fespace as an observer of mesh
@@ -48,19 +50,24 @@ classdef FeSpace < handle
obj.mesh = mesh;
obj.finiteElement = fe;
- % store indices of dirichlet boundary
- if isa(opt.dirichlet, 'char')
- if strcmp(opt.dirichlet, 'all')
- obj.dirichlet = 1:mesh.nBoundaries;
- else
- error('FeSpace:invalidDirichletBoundary', 'Choice of dirichlet boundary not valid!')
- end
- elseif isa(opt.dirichlet, 'double')
- if any(opt.dirichlet < 1 | opt.dirichlet > mesh.nBoundaries, 'all')
- error('FeSpace:invalidDirichletBoundary', 'Choice of dirichlet boundary not valid!')
- else
- obj.dirichlet = unique(opt.dirichlet);
+ % store indices of boundaries and do sanity checks
+ obj.bnd = bnd;
+ if all(structfun(@isempty, bnd))
+ obj.bnd.dirichlet = ':';
+ end
+
+ bndNames = {'dirichlet', 'neumann', 'robin'};
+ referenceCount = zeros(1, mesh.nBoundaries);
+ for k=1:numel(bndNames)
+ idx = obj.bnd.(bndNames{k});
+ if ~ischar(idx)
+ mustBeInRange(idx, 1, mesh.nBoundaries)
end
+ referenceCount(idx) = referenceCount(idx) + 1;
+ end
+ if any(referenceCount ~= 1)
+ error('FeSpace:invalidBoundary', ...
+ 'Every boundary part must be connected to exactly one boundary.')
end
end
@@ -77,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/assembleFreeDofs.m b/lib/spaces/@FeSpace/assembleFreeDofs.m
index 6637ed7bb9323ce2cb37344d48870ce2ea722ce3..864617d392e2d22f64ed67db5adecd087d5008c3 100644
--- a/lib/spaces/@FeSpace/assembleFreeDofs.m
+++ b/lib/spaces/@FeSpace/assembleFreeDofs.m
@@ -8,7 +8,7 @@ dofs = getDofs(obj);
mesh = obj.mesh;
allDofs = true(dofs.nDofs, 1);
-dirichletEdges = getCombinedBndEdges(mesh, obj.dirichlet);
+dirichletEdges = getCombinedBndEdges(mesh, obj.bnd.dirichlet);
allDofs(dofs.edge2Dofs(:,dirichletEdges)) = false;
freeDofs = find(allDofs);
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/tests/TestAssembly.m b/tests/TestAssembly.m
index d34353a59e9909cddb8c281af13c37ff083ac052..9a9a40b1024cad74827361d70f521fe06654d01e 100644
--- a/tests/TestAssembly.m
+++ b/tests/TestAssembly.m
@@ -2,6 +2,7 @@ classdef TestAssembly < matlab.unittest.TestCase
properties
mesh
+ fes
blf
lf
end
@@ -10,82 +11,81 @@ methods (TestMethodSetup)
function initializeForms(testCase)
% expected values below are computed by hand for unit triangle!
testCase.mesh = Mesh.loadFromGeometry('unittriangle');
- fes = FeSpace(testCase.mesh, LowestOrderH1Fe());
- testCase.blf = BilinearForm(fes);
- testCase.lf = LinearForm(fes);
+ testCase.fes = FeSpace(testCase.mesh, LowestOrderH1Fe());
+ testCase.blf = BilinearForm();
+ testCase.lf = LinearForm();
end
end
methods (Test)
function assembleWithoutCoefficientThrowsError(testCase)
testCase.blf.a = [];
- testCase.verifyError(@() assemble(testCase.blf), 'BilinearForm:NoCoefficients')
- testCase.verifyError(@() assemble(testCase.lf), 'LinearForm:NoCoefficients')
+ testCase.verifyError(@() assemble(testCase.blf, testCase.fes), 'BilinearForm:NoCoefficients')
+ testCase.verifyError(@() assemble(testCase.lf, testCase.fes), 'LinearForm:NoCoefficients')
end
function validateStiffnessScalar(testCase)
testCase.blf.a = Constant(testCase.mesh, 1);
expected = [2,-1,-1;-1,1,0;-1,0,1]/2;
- testCase.verifyEqual(full(assemble(testCase.blf)), expected, 'RelTol', eps);
+ testCase.verifyEqual(full(assemble(testCase.blf, testCase.fes)), expected, 'RelTol', eps);
end
function validateStiffnessMatrix(testCase)
testCase.blf.a = Constant(testCase.mesh, [1;0;1;2]);
expected = [4,-1,-3;-2,1,1;-2,0,2]/2;
- testCase.verifyEqual(full(assemble(testCase.blf)), expected, 'RelTol', eps);
+ testCase.verifyEqual(full(assemble(testCase.blf, testCase.fes)), expected, 'RelTol', eps);
end
function validateConvection(testCase)
testCase.blf.a = [];
testCase.blf.b = Constant(testCase.mesh, [1;1]);
expected = [-2,1,1;-2,1,1;-2,1,1]/6;
- testCase.verifyEqual(full(assemble(testCase.blf)), expected, 'RelTol', eps);
+ testCase.verifyEqual(full(assemble(testCase.blf, testCase.fes)), expected, 'RelTol', eps);
end
function validateReactionWithHigherOrderQuadrature(testCase)
testCase.blf.c = MeshFunction(testCase.mesh, @(x) 120*sum(x, Dim.Vector));
expected = [4,3,3;3,8,4;3,4,8];
- testCase.verifyNotEqual(full(assemble(testCase.blf)), expected);
+ testCase.verifyNotEqual(full(assemble(testCase.blf, testCase.fes)), expected);
testCase.blf.qrc = QuadratureRule.ofOrder(2);
- testCase.verifyNotEqual(full(assemble(testCase.blf)), expected);
+ testCase.verifyNotEqual(full(assemble(testCase.blf, testCase.fes)), expected);
testCase.blf.qrc = QuadratureRule.ofOrder(3);
- testCase.verifyEqual(full(assemble(testCase.blf)), expected, 'RelTol', eps*10);
+ testCase.verifyEqual(full(assemble(testCase.blf, testCase.fes)), expected, 'RelTol', eps*10);
end
function validateL2rhsScalar(testCase)
testCase.lf.f = Constant(testCase.mesh, 1);
expected = [1;1;1]/6;
- testCase.verifyEqual(assemble(testCase.lf), expected, 'RelTol', eps);
+ testCase.verifyEqual(assemble(testCase.lf, testCase.fes), expected, 'RelTol', eps);
end
function validateH1rhsMeshFunctionWithHigherOrderQuadrature(testCase)
testCase.lf.fvec = MeshFunction(testCase.mesh, ...
@(x) [120;0].*(prod(x, Dim.Vector).*(1-sum(x, Dim.Vector))));
expected = [-1;1;0];
- testCase.verifyNotEqual(assemble(testCase.lf), expected);
+ testCase.verifyNotEqual(assemble(testCase.lf, testCase.fes), expected);
testCase.lf.qrfvec = QuadratureRule.ofOrder(2);
- testCase.verifyNotEqual(assemble(testCase.lf), expected);
+ testCase.verifyNotEqual(assemble(testCase.lf, testCase.fes), expected);
testCase.lf.qrfvec = QuadratureRule.ofOrder(3);
- testCase.verifyEqual(assemble(testCase.lf), expected, 'RelTol', eps*10);
+ testCase.verifyEqual(assemble(testCase.lf, testCase.fes), expected, 'RelTol', eps*10);
end
function validateNeumann(testCase)
- testCase.lf.neumann = Constant(testCase.mesh, pi);
- testCase.lf.bndNeumann = [1,2];
+ neumannFes = FeSpace(testCase.mesh, LowestOrderH1Fe(), 'neumann', ':');
+ nlf = LinearForm();
+ nlf.neumann = Constant(testCase.mesh, pi);
expected = [2;1+sqrt(2);1+sqrt(2)]*(pi/2);
- testCase.verifyEqual(assemble(testCase.lf), expected, 'RelTol', eps*10);
+ testCase.verifyEqual(assemble(nlf, neumannFes), expected, 'RelTol', eps*10);
end
function validateRobin(testCase)
- testCase.lf.neumann = Constant(testCase.mesh, pi);
- testCase.lf.bndNeumann = ':';
- expected = [2;1+sqrt(2);1+sqrt(2)]*(pi/2);
- testCase.verifyEqual(assemble(testCase.lf), expected, 'RelTol', eps*10);
- testCase.blf.robin = Constant(testCase.mesh, pi);
- testCase.blf.bndRobin = 1;
- testCase.blf.qrRobin = QuadratureRule.ofOrder(2, '1D');
+ robinFes = FeSpace(testCase.mesh, LowestOrderH1Fe(), 'dirichlet', 2, 'robin', 1);
+ rlf = LinearForm();
+ rblf = BilinearForm();
+ rblf.robin = Constant(testCase.mesh, pi);
+ rblf.qrRobin = QuadratureRule.ofOrder(2, '1D');
expected = [4,1,1;1,2,0;1,0,2]*(pi/6);
- testCase.verifyEqual(full(assemble(testCase.blf)), expected, 'RelTol', eps*10);
+ testCase.verifyEqual(full(assemble(rblf, robinFes)), expected, 'RelTol', eps*10);
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 48cef84fe9bec01420268c91a1c4ba8c9a563009..fa8f58e1c70e52025fb4f913b77c64600ab4f40a 100644
--- a/tests/TestFeSpace.m
+++ b/tests/TestFeSpace.m
@@ -5,8 +5,7 @@ properties
end
properties (TestParameter)
- dirichlet = {1, 2, [1,2], [1,2,1,2,1,1], 'all'};
- expectedDirichlet = {1, 2, [1,2], [1,2], [1,2]};
+ dirichlet = {[1,2], [2,1], ':'};
fe = struct('H1p1', LowestOrderH1Fe(), 'L2p1', LowestOrderL2Fe(), ...
'H1p3', HigherOrderH1Fe(4), 'L2p7', HigherOrderL2Fe(7));
nDof = {3, 1, 15, 36};
@@ -23,27 +22,41 @@ methods (Test, ParameterCombination='sequential')
function implicitDirichletTakesAllBoundaries(testCase)
dummyFe = LowestOrderH1Fe();
fes = FeSpace(testCase.dummyMesh, dummyFe);
- testCase.verifyEqual(fes.dirichlet, [1,2])
+ boundaries = 1:fes.mesh.nBoundaries;
+ testCase.verifyEqual(boundaries(fes.bnd.dirichlet)', [1,2])
end
- function allOptionsForDirichletGiveCorrectBoundaries(testCase, dirichlet, expectedDirichlet)
+ function allOptionsForDirichletGiveCorrectBoundaries(testCase, dirichlet)
dummyFe = LowestOrderH1Fe();
fes = FeSpace(testCase.dummyMesh, dummyFe, 'dirichlet', dirichlet);
- testCase.verifyEqual(fes.dirichlet, expectedDirichlet);
+ testCase.verifyEqual(fes.bnd.dirichlet, dirichlet);
end
- function wrongDirichletDataThrowsError(testCase)
+ function invalidBoundaryDataThrowsError(testCase)
dummyFe = LowestOrderH1Fe();
testCase.verifyError( ...
- @() FeSpace(testCase.dummyMesh, dummyFe, 'dirichlet', [1,3,1,2,1,1]), ...
- 'FeSpace:invalidDirichletBoundary')
+ @() FeSpace(testCase.dummyMesh, dummyFe, 'dirichlet', [1,2,3]), ...
+ 'MATLAB:validators:mustBeInRange')
testCase.verifyError( ...
- @() FeSpace(testCase.dummyMesh, dummyFe, 'dirichlet', 'whole'), ...
- 'FeSpace:invalidDirichletBoundary')
+ @() FeSpace(testCase.dummyMesh, dummyFe, 'neumann', 'whole'), ...
+ 'mustBeIndexVector:notAnIndexVector')
+ end
+
+ function boundaryPartsMakeUpWholeBoundary(testCase)
+ dummyFe = LowestOrderH1Fe();
+ testCase.verifyError( ...
+ @() FeSpace(testCase.dummyMesh, dummyFe, 'robin', 1), ...
+ 'FeSpace:invalidBoundary')
+ testCase.verifyError( ...
+ @() FeSpace(testCase.dummyMesh, dummyFe, 'neumann', 1, 'robin', 1), ...
+ 'FeSpace:invalidBoundary')
+ testCase.verifyWarningFree(...
+ @() FeSpace(testCase.dummyMesh, dummyFe, ...
+ 'dirichlet', [], 'neumann', 1, 'robin', 2))
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..11964911c3bc50586be4a40a1b716d06cf1dc175
--- /dev/null
+++ b/tests/TestIterativeSolver.m
@@ -0,0 +1,92 @@
+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.c = Constant(mesh, 1);
+ 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/TestMesh.m b/tests/TestMesh.m
index ea0e2972e4e702dbab7c039cd92a5c653d5daa0e..225ff5d692fa02b689f3e9ad7c89c225a9a5c2ee 100644
--- a/tests/TestMesh.m
+++ b/tests/TestMesh.m
@@ -11,6 +11,17 @@ methods (Test)
testCase.verifyClass(Mesh.loadFromGeometry(folders{i}), 'Mesh')
end
end
+
+ function cloneMethodMakesClone(testCase)
+ mesh = Mesh.loadFromGeometry('Lshape');
+ clonedMesh = clone(mesh);
+ testCase.verifyEqual(mesh.coordinates, clonedMesh.coordinates);
+ testCase.verifyEqual(mesh.elements, clonedMesh.elements);
+ testCase.verifyEqual(mesh.edges, clonedMesh.edges);
+ testCase.verifyEqual(mesh.element2edges, clonedMesh.element2edges);
+ testCase.verifyEqual(mesh.flipEdges, clonedMesh.flipEdges);
+ testCase.verifyEqual(mesh.boundaries, clonedMesh.boundaries);
+ end
function LShapeDetailsAreCorrect(testCase)
% this is tailored to the Lshape to ensure correctness of arrays
diff --git a/tests/TestPatchwiseMatrix.m b/tests/TestPatchwiseMatrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..206b6a6351075c53ed9b1797c5e783ed868cbaac
--- /dev/null
+++ b/tests/TestPatchwiseMatrix.m
@@ -0,0 +1,68 @@
+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.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