diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
index f353744f765d7017efa7c02823f5d3f8c96a3b8a..9e6211ebd724ca7f80ce5d711e49ab5a679322c5 100644
--- a/examples/algebraicSolver.m
+++ b/examples/algebraicSolver.m
@@ -30,7 +30,9 @@ for p = 1:pmax
% choose the iterative solver of choice: multigrid (with the variants
% lowOrderVcycle and highOrderVcycle), pcg (with jacobi and additive
% Schwarz preconditioner), and the cg solver
- [solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
+ %[solver, 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;
diff --git a/experiments/p_laplace.m b/experiments/p_laplace.m
deleted file mode 100644
index 16e81435101ff679cd305a62701dcc28128c05f7..0000000000000000000000000000000000000000
--- a/experiments/p_laplace.m
+++ /dev/null
@@ -1,170 +0,0 @@
-%% paramters
-nElemMax = 1e5;
-theta = 0.5;
-lambda = 0.1;
-deltaZ = 0.1;
-linearization = 'kacanov';
-geometry = 'Lshape';
-p = 3;
-
-%% initialisation
-[nElem, eta, E] = deal(zeros(1,1000));
-
-%% setup geometry & spaces
-mesh = Mesh.loadFromGeometry(geometry);
-mesh.refineUniform(2);
-fes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', ':');
-u = FeFunction(fes);
-
-%% set problem data given in [Heid, Praetorius, Wihler; 2021]
-Du = Gradient(u);
-g = Constant(mesh, 1);
-[blf, lf] = setProblemData(mesh, fes, Du, g, linearization, p);
-eDensity = CompositeFunction(@(t) muIntegral(vectorProduct(t, t),p), Du);
-u.setFreeData(rand(1,length(getFreeDofs(u.fes))));
-
-%% nested iteration
-P = LoMeshProlongation(fes);
-
-%% adaptive algorithm
-i = 1;
-Eold = 0;
-qr = QuadratureRule.ofOrder(1); % tweak?
-while 1
- % the volume residual for lowest order is independent of the approximate
- % solution, so it can be computed once
- vol = estimateVolume(g);
-
- %% inner iteration to solve the nonlinear problem
- nIterations = 1;
- while true
- %% solve linearized problem
- A = assemble(blf, fes);
- F = assemble(lf, fes);
- updateDataU(u, deltaZ, A, F, linearization);
-
- %% compute remainder of error estimator and energy
- edge = estimateEdge(u, mesh, p);
- eta2 = combineEstimators(vol, edge, mesh);
- eta(i) = sqrt(sum(eta2));
- nElem(i) = mesh.nElements;
-
- if i > 1, Eold = E(i-1); end
- E(i) = sum(integrateElement(eDensity, qr)) + u.data * F;
-
- %% stopping criterion of the iterative solver
- if sqrt(E(i) - Eold) <= lambda * eta(i) % norm(e) <= fehler
- break
- end
- nIterations = nIterations + 1;
- i = i+1;
- end
-
- %% stoping criterion for the AFEM algorithm
- printLogMessage('number of elements: %d, iterations: %d, estimator: %.2e', ...
- nElem(i), nIterations, eta(i));
- if nElem(i) > nElemMax
- break
- end
-
- %% refine mesh
- marked = markDoerflerSorting(eta2, theta);
- mesh.refineLocally(marked, 'NVB');
- u.setData(prolongate(P, u));
- i = i+1;
-end
-
-%% plot convergence rates
-plotData(nElem, 'number of elements', eta, 'error estimator', linearization,p)
-
-%% helper function for plotting
-function plotData(xdata, xlab, ydata, ylab, name,p)
- figure()
- for k = 1:size(xdata,1)
- idx = find(xdata(k,:) > 0);
- loglog(xdata(k,idx), ydata(k,idx), '-o', 'LineWidth', 2, 'DisplayName', strcat(name,', p=', num2str(p)));
- hold on
- end
- x = xdata(1,idx) / xdata(1,1);
- loglog(xdata(1,1)*x, ydata(1,1)*x.^(-1/2), '--', 'LineWidth', 2, 'Color', 'k', 'DisplayName', '\alpha = 1/2')
- legend
- xlabel(xlab)
- ylabel(ylab)
- title([ylab, ' over ', xlab])
-end
-
-%% local function for problem data
-function [blf, lf] = setProblemData(mesh, fes, Du, g, linearization, p)
- blf = BilinearForm();
- lf = LinearForm();
-
- switch linearization
- case "zarantonello"
- blf.a = Constant(mesh, 1);
- lf.f = g;
- lf.fvec = CompositeFunction(@(t) -mu(vectorProduct(t, t),p) .* t, Du);
- case "kacanov"
- blf.a = CompositeFunction(@(t) mu(vectorProduct(t, t),p), Du);
- lf.f = g;
- lf.fvec = [];
- case "newton"
- blf.a = CompositeFunction(@(t) mu(vectorProduct(t, t),p) .* [1;0;0;1] ...
- + 2*muDerivative(vectorProduct(t, t),p).*vectorProduct(t, t, [2,1], [2,1]'), Du);
- lf.f = g;
- lf.fvec = CompositeFunction(@(t) -mu(vectorProduct(t, t),p) .* t, Du);
- end
-end
-
-function updateDataU(u, deltaZ, A, F, linearization)
- freeDofs = getFreeDofs(u.fes);
- v = FeFunction(u.fes);
- switch linearization
- case "zarantonello"
- v.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
- u.setData(u.data + deltaZ*v.data);
- case "kacanov"
- %u.setData(A \ F);
- u.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
- case "newton"
- v.setFreeData(A(freeDofs,freeDofs) \ F(freeDofs));
- u.setData(u.data + v.data);
- end
-end
-
-%% local function for residual a posteriori error estimation
-% \eta(T)^2 = h_T^2 * || g ||_{L^2(T)}^2
-% + h_T * || [[ mu(|Du|^2)*Du * n]] ||_{L^2(E) \cap \Omega}^2
-function vol = estimateVolume(g)
- qr = QuadratureRule.ofOrder(1);
- vol = integrateElement(g, qr);
-end
-
-function edge = estimateEdge(u, mesh, p)
- qr = QuadratureRule.ofOrder(1, '1D');
- dirichlet = vertcat(mesh.boundaries{:});
-
- f = CompositeFunction(@(t) mu(vectorProduct(t, t),p).*t, Gradient(u));
- edge = integrateNormalJump(f, qr, ...
- @(j) zeros(size(j)), {}, dirichlet, ... % no jump on dirichlet edges
- @(j) j.^2, {}, ':'); % normal jump on inner edges
-
-end
-
-function indicators = combineEstimators(vol, edge, mesh)
- trafo = getAffineTransformation(mesh);
- hT = sqrt(trafo.area);
- indicators = hT.^2.*vol + hT.*sum(edge(mesh.element2edges), Dim.Vector);
-end
-
-%% nonlinearities
-function y = mu(t,p)
- y = t.^(p/2 - 1);
-end
-
-function y = muIntegral(t,p)
- y = p/4 * t.^(p/2);
-end
-
-function y = muDerivative(t,p)
- y = (p/2 - 1) * t.^(p/2 - 2);
-end
diff --git a/runAdditiveSchwarzPreconditioner.m b/runAdditiveSchwarzPreconditioner.m
deleted file mode 100644
index a0dbc376bbeba500ba22ff0476058fae1069bdce..0000000000000000000000000000000000000000
--- a/runAdditiveSchwarzPreconditioner.m
+++ /dev/null
@@ -1,123 +0,0 @@
-function leveldata = runAdditiveSchwarzPreconditioner(pMax)
-
- % Proceed input
- if nargin < 1
- pMax = 1e1;
- end
-
- % Number of refinement levels
- nLevels = 2;
-
- % Initialize LevelData
- leveldata = LevelData('results');
- leveldata.problem = 'Poisson';
- leveldata.domain = 'UnitSquare';
- leveldata.method = 'Sp_PCG_AdditiveSchwarz';
-
- % Run experiment for various p
- for p = 1:pMax
- % Load mesh
- mesh = Mesh.loadFromGeometry('unitsquare');
- %mesh.refineUniform();
-
- % Create FE space
- fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', ':');
- % u = FeFunction(fes);
-
- % Setup variational problem
- blf = BilinearForm();
- lf = LinearForm();
- blf.a = Constant(mesh, 1);
- lf.f = Constant(mesh, 1);
-
- % Choose iterative solver
- solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
- % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
- % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
- % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
-
- for k = 1:nLevels-1
- % Assemble matrices on intermediate meshes
- A = assemble(blf, fes);
- F = assemble(lf, fes);
- freeDofs = getFreeDofs(fes);
- solver.setupSystemMatrix(A(freeDofs,freeDofs));
- solver.setupRhs(F(freeDofs)); % initializes x0 as well
- % Mesh refinement
- mesh.refineUniform();
- end
-
- % Assemble matrices on fine mesh
- A = assemble(blf, fes);
- if ~issymmetric(A)
- %warning(['symmetrizing, norm(asym(A)) = ', num2str(norm(A - A', 1))])
- A = (A + A') / 2;
- end
- F = assemble(lf, fes);
- freeDofs = getFreeDofs(fes);
- solver.setupSystemMatrix(A(freeDofs,freeDofs));
- solver.setupRhs(F(freeDofs)); % initializes x0 as well
-
- % A = A(freeDofs,freeDofs);
-
- % figure(2);
- % spy(A);
- % title('Matrix A');
-
- % Compute C
- % C = zeros(length(freeDofs));
- % for j = 1:length(freeDofs)
- % ej = zeros(length(freeDofs), 1);
- % ej(j) = 1;
- % C(:,j) = solver.preconditionAction(ej);
- % end
- % if ~issymmetric(C)
- % % warning(['symmetrizing, norm(asym(C)) = ', num2str(norm(C - C', 1))])
- % C = (C + C') / 2;
- % end
-
- % figure(3);
- % CC = C;
- % CC(abs(C) < 1e-10) = 0;
- % spy(CC);
- % title('Matrix C');
- %
- % asymC = (C - C');
-
- % figure(4);
- % largeAsymC = asymC;
- % largeAsymC(abs(largeAsymC) < max(abs(largeAsymC(:))) / 10) = 0;
- % spy(largeAsymC);
- % title('Large entries of asym(C)');
-
- % figure(5);
- % asymC(abs(asymC) < 1e-10) = 0;
- % spy(asymC);
- % title('Nonzero entries of asym(C)');
-
- % Compute conditional number
- % CA = C*A;
- % lambda = eig(full(CA));
- % condition = max(abs(lambda)) / min(abs(lambda));
- % lambdaMin = eigs(@(x)applyCA(x, A, solver), length(freeDofs), 1, 'smallestreal');
- % lambdaMax = eigs(@(x)applyCA(x, A, solver), length(freeDofs), 1, 'largestreal');
- % lambdaMin = eigs(CA, 1, 'smallestabs');
- % lambdaMax = eigs(CA, 1, 'largestabs');
- % condition = lambdaMax / lambdaMin;
-
- % Print result to commandline
- % leveldata.append('p', uint32(p), 'condition', condition);
- % leveldata.printLevel();
- end
-
- % Plot condition
- figure();
- leveldata.plot('p');
-end
-
-
-function CAx = applyCA(x, A, solver)
- % CAx = solver.preconditionAction(A * x);
- % CAx = (A * x);
- CAx = solver.preconditionAction(x);
-end
diff --git a/runIterativeSolver.m b/runIterativeSolver.m
deleted file mode 100644
index a14ceae87e6bf7a7b550c59b7b2b0b33b338c0ec..0000000000000000000000000000000000000000
--- a/runIterativeSolver.m
+++ /dev/null
@@ -1,144 +0,0 @@
-function leveldata = runIterativeSolver(maxNiter)
-
- % Proceed input
- if nargin < 1
- maxNiter = 1e2;
- end
-
- % Polynomial degree
- p = 4;
-
- % Number of refinement levels
- nLevels = 6;
-
- % Load mesh
- mesh = Mesh.loadFromGeometry('unitsquare');
- mesh.refineUniform();
-
- % Create FE space
- fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', ':');
- u = FeFunction(fes);
- ustar = FeFunction(fes);
-
- % Setup variational problem
- blf = BilinearForm();
- lf = LinearForm();
- blf.a = Constant(mesh, 1);
- lf.f = Constant(mesh, 1);
-
- % Choose iterative solver
- solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
- % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
- % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
- % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
-
- % Initialize LevelData
- leveldata = LevelData('results');
- leveldata.problem = 'Poisson';
- leveldata.domain = 'UnitSquare';
- leveldata.method = 'Sp_PCG_AdditiveSchwarz';
-
- for k = 1:nLevels-1
- % Output
- fprintf('Assemble data for level = %d\n', k);
- % Assemble matrices on intermediate meshes
- A = assemble(blf, fes);
- F = assemble(lf, fes);
- freeDofs = getFreeDofs(fes);
- solver.setupSystemMatrix(A(freeDofs,freeDofs));
- solver.setupRhs(F(freeDofs)); % initializes x0 as well
- % Mesh refinement
- mesh.refineUniform();
- end
-
- fprintf('Assemble data for level = %d\n', nLevels);
-
- % Assemble matrices on fine mesh
- A = assemble(blf, fes);
- F = assemble(lf, fes);
- freeDofs = getFreeDofs(fes);
- solver.setupSystemMatrix(A(freeDofs,freeDofs));
- solver.setupRhs(F(freeDofs)); % initializes x0 as well
-
- % Compute exact solution
- fprintf('Do exact solve\n');
- xstar = A(freeDofs,freeDofs) \ F(freeDofs);
-
- % Refinement loop
- nIter = 1;
- while true
- ndof = getDofs(fes).nDofs;
-
- % Compute energy norm
- errorU = sqrt((xstar - solver.x)' * A(freeDofs,freeDofs) * (xstar - solver.x));
-
- % Compute contraction factor for energy norm
- if nIter == 1
- contraction = Inf;
- else
- contraction = errorU / errorUPrevious;
- end
-
- % Store variables
- leveldata.append('ndof', uint32(ndof), ...
- 'errorU', errorU);
- leveldata.setAbsolute(nIter, 'nIter', uint32(nIter), ...
- 'contraction', contraction);
-
- % Print level information
- leveldata.printLevel();
-
- % Plot solution
- % if p == 2
- % x = zeros(mesh.nCoordinates + mesh.nEdges, 1);
- % x(freeDofs) = solver.x;
- % figure(1);
- % plotS2(mesh.clone(), x);
- %
- % figure(2);
- % x(freeDofs) = xstar;
- % plotS2(mesh.clone(), x);
- % title("Exact algebraic solution")
- % else
- % figure(1);
- % u.setFreeData(solver.x);
- % plot(u);
- %
- % figure(2);
- % ustar.setFreeData(xstar);
- % plot(ustar);
- % end
- % pause(2)
-
- % Break condition
- % if ndof > maxNdof
- % break;
- % end
-
- % Mesh refinement
- % mesh.refineUniform();
-
- if nIter > maxNiter
- break;
- end
-
- % One iterative solver step
- solver.step();
-
- % Update iteration counter
- nIter = nIter + 1;
- errorUPrevious = errorU;
- end
-
- % Plot solution
- u.setFreeData(solver.x);
- plot(u);
-
- % Plot error
- figure();
- leveldata.plot('nIter');
-
- % Plot contraction
- figure();
- leveldata.plotAbsolute('nIter');
-end
diff --git a/test.m b/test.m
deleted file mode 100644
index 505e85f64341a4cb0f5b3508fa0e7e2d6eb4d257..0000000000000000000000000000000000000000
--- a/test.m
+++ /dev/null
@@ -1,61 +0,0 @@
-mesh = Mesh.loadFromGeometry('unitsquare');
-%mesh.refineUniform();
-
-fes1 = FeSpace(mesh, HigherOrderH1Fe(1), 'dirichlet', ':');
-fes2 = FeSpace(mesh, HigherOrderH1Fe(15), 'dirichlet', ':');
-
-testMatrix = ones(getDofs(fes2).nDofs, 1);
-inclusionMatrix = MeshProlongation(fes1, fes2);
-inclusionMatrix = full(inclusionMatrix);
-Vertices = 1:getDofs(fes1).nDofs;
-inclusionMatrix2 = interpolateData(eye(length(Vertices)), fes1, fes2)';
-diff = inclusionMatrix - inclusionMatrix2;
-
-function inclusionMatrix = MeshProlongation(fromFes, toFes)
- %Create the prolongation matrix on the unit triangle
- unittriangle = Mesh.loadFromGeometry('unittriangle');
- fromFesUnitTriangle = FeSpace(unittriangle, ...
- HigherOrderH1Fe(fromFes.finiteElement.order));
- toFesUnitTriangle = FeSpace(unittriangle, ...
- HigherOrderH1Fe(toFes.finiteElement.order));
- Vertices = 1:getDofs(fromFesUnitTriangle).nDofs;
- unitTriangleInclusionMatrix = ...
- permute(interpolateData(eye(length(Vertices)),...
- fromFesUnitTriangle, toFesUnitTriangle), [2 1]);
- % Get local to global map in unit triangle
-
- unitTriangleInclusionMatrix = unitTriangleInclusionMatrix( ...
- getDofs(fromFesUnitTriangle).element2Dofs, ...
- getDofs(toFesUnitTriangle).element2Dofs);
- % Get dofs of the finite element spaces on the meshes
- fromFesDofs = getDofs(fromFes);
- toFesDofs = getDofs(toFes);
-
- % Copy the matrix from the unit triangle for each element
- mat = repmat(unitTriangleInclusionMatrix, ...
- [1, 1, fromFes.mesh.nElements]);
-
- % Create index sets for the accumarray
- temp = fromFesDofs.element2Dofs(:, :);
- temp2 = toFesDofs.element2Dofs(:, :);
- I = repmat(permute(temp, [1 3 2]), [1, size(temp2, 1), 1]);
- J = repmat(permute(temp2, [3 1 2]), [size(temp, 1), 1, 1]);
- ind = [I(:), J(:)];
-
- inclusionMatrix = ...
- accumarray(ind, mat(:), [], @mean, [], true);
-end
-
-function interpolatedData = interpolateData(data, fromFes, toFes)
- Dofs = 1:getDofs(toFes).nDofs;
- nComponents = size(data, 2);
- interpolatedData = zeros(numel(Dofs), nComponents);
- feFunctionWrapper = FeFunction(fromFes);
- for k = 1:nComponents
- feFunctionWrapper.setData(data(:,k));
- wholeData = nodalInterpolation(feFunctionWrapper, toFes);
- interpolatedData(:,k) = wholeData(Dofs);
- end
-end
-
-
diff --git a/tmp/plotS2.m b/tmp/plotS2.m
deleted file mode 100644
index 92ee44046520096223c78bb8688872588c419297..0000000000000000000000000000000000000000
--- a/tmp/plotS2.m
+++ /dev/null
@@ -1,49 +0,0 @@
-function ax = plotS2(mesh, x)
-%%PLOTS2 plots a piecewise P2 and globally continuous function (Lagrange
-%finite element function) given by the coefficient vector x on the
-%triangulation meshdata
-% ax = PLOTS2(meshdata, x)
-
-% Copyright 2023 Philipp Bringmann
-%
-% This program is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% This program is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-%
-
-
- % Create axes object
- ax = gca;
-
- % Uniform red-refinement of triangulation
- mesh.refineUniform();
-
- % S1 plot of nodal interpolation on red-refinement
- if size(x, 2) == 1
- % Plot of scalar function
- pat = trisurf(mesh.elements', mesh.coordinates(1,:)',...
- mesh.coordinates(2,:)', x);
- if(mesh.nElements > 2000)
- set(pat, 'EdgeColor', 'none');
- end
- elseif size(x, 2) == 2
- % Plot of vector field
- set(ax, 'XLim', 2.1*[min(mesh.coordinates(1,:)), max(mesh.coordinates(1,:))],...
- 'YLim', 1.1*[min(mesh.coordinates(2,:)), max(mesh.coordinates(2,:))]);
- quiver(mesh.coordinates(1,:), mesh.coordinates(2,:), x(:,1), x(:,2));
- else
- error('Invalid size of coefficient vector');
- end
-
- % Draw title
- title(ax, {'S2 function plot'; [num2str(mesh.nEdges + mesh.nCoordinates), ' dofs']});
-end