diff --git a/lib/solvers/LocalMgLowOrderVcycle.m b/lib/solvers/LocalMgLowOrderVcycle.m
deleted file mode 100644
index f8920c03d1057f10ed21b3f411112f4e4c83abd5..0000000000000000000000000000000000000000
--- a/lib/solvers/LocalMgLowOrderVcycle.m
+++ /dev/null
@@ -1,179 +0,0 @@
-% OptimalLocalMGSolver (subclass of MGSolver) Solves linear equations
-% iteratively using high-order p-optimal multilevel local smoothing (additive
-% Schwarz/block Jacobi).
-% One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
-%
-% With the solver comes an
-
-classdef LocalMgLowOrderVcycle < MGSolver
- %% properties
- properties (GetAccess=public,SetAccess=protected)
- nLevels
- end
-
- properties (Access=protected)
- blf
- hoFes
- loFes
- p1Matrix
- p1Smoother
- P
- intergridMatrix
- inclusionMatrix
- freeVertices
- freeVerticesOld
- changedPatches
- patchwiseA
- end
-
- %% event data
- properties (Access=private)
- listenerHandle
- end
-
- %% methods
- methods (Access=public)
- function obj = LocalMgLowOrderVcycle(fes, blf)
- arguments
- fes FeSpace
- blf BilinearForm
- end
- obj = obj@MGSolver();
-
- assert(fes.finiteElement.order > 1, ...
- 'LocalMgLowOrderVcycle only works for higher order finite elements.')
- assert(isempty(blf.b), ...
- 'Multigrid solvers only tested for symmetric problems.')
-
- obj.nLevels = 0;
-
- mesh = fes.mesh;
- obj.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', ':');
- obj.hoFes = fes;
- obj.P = Prolongation.chooseFor(obj.loFes);
- obj.blf = blf;
-
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
- end
-
- function setupSystemMatrix(obj, A)
- obj.nLevels = obj.nLevels + 1;
-
- L = obj.nLevels;
- obj.freeVerticesOld = obj.freeVertices;
- obj.freeVertices = getFreeDofs(obj.loFes);
- freeVerticesHigher = getFreeDofs(obj.hoFes);
-
- obj.p1Matrix{L} = assemble(obj.blf, obj.loFes);
- obj.p1Matrix{L} = obj.p1Matrix{L}(obj.freeVertices, obj.freeVertices);
- obj.p1Smoother{L} = full(diag(obj.p1Matrix{L})).^(-1);
-
- if L >= 2
- obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
- obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
- obj.patchwiseA = assemblePatchwise(obj.blf, obj.hoFes);
- obj.inclusionMatrix = SpaceProlongation(obj.loFes, obj.hoFes) ...
- .matrix(getFreeDofs(obj.loFes), getFreeDofs(obj.hoFes));
- end
- setupSystemMatrix@MGSolver(obj, A);
- end
-
- function setupRhs(obj, b, varargin)
- setupRhs@MGSolver(obj, b, varargin{:});
- end
-
- % Geometric MultiGrid
- function [Cx, algError2] = Vcycle(obj, res)
- assert(isequal(size(res, 1), size(obj.A, 1)), ...
- 'Setup for multilevel iteration not complete!')
-
- % if there is only one coarse level: exact solve
- if obj.nLevels == 1
- Cx = obj.A \ res;
- algError2 = scalarProduct(Cx, obj.A*Cx);
- return
- end
-
- L = obj.nLevels;
- algError2 = zeros(1, size(res, 2)); % built-in estimator of the algebraic error
-
- % descending cascade in P1: no smoothing
- residual{L} = obj.inclusionMatrix * res; %interpolateFreeData(res, obj.hoFes, obj.loFes);
- for k = L:-1:2
- residual{k-1} = obj.intergridMatrix{k}'*residual{k};
- end
-
- % exact solve on coarsest level to compute accumulative lifting
- % of the residual (sigma)
- sigma = obj.p1Matrix{1} \ residual{1};
- algError2 = algError2 + scalarProduct(sigma, obj.p1Matrix{1}*sigma);
-
- % ascending cascade in P1 AND local smoothing
- for k = 2:(L-1)
- sigma = obj.intergridMatrix{k}*sigma;
- uptres = residual{k} - obj.p1Matrix{k}*sigma; %updated residual
- rho = p1LocalSmoothing(obj, k, uptres);
- [aeUpd, sUpd] = computeOptimalUpdate(obj.p1Matrix{k}, uptres, rho);
- sigma = sigma + sUpd;
- algError2 = algError2 + aeUpd;
- end
-
- % smoothing on finest level dependent on p
- sigma = obj.intergridMatrix{L}*sigma;
- sigma = obj.inclusionMatrix' * sigma; %interpolateFreeData(sigma, obj.loFes, obj.hoFes);
- uptres = res - obj.A*sigma;
- rho = hoGlobalSmoothing(obj, uptres);
- [aeUpd, sUpd] = computeOptimalUpdate(obj.A, uptres, rho);
- sigma = sigma + sUpd;
- algError2 = algError2 + aeUpd;
-
- Cx = sigma;
- 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
-
- function rho = p1LocalSmoothing(obj, k, res)
- idx = obj.changedPatches{k};
- rho = zeros(size(res));
- rho(idx,:) = obj.p1Smoother{k}(idx).*res(idx,:);
- end
-
- function rho = hoGlobalSmoothing(obj, res)
- rho = obj.patchwiseA \ res;
- end
- end
-end
-
-%% auxiliary functions
-% error correction with optimal stepsize
-function [etaUpdate, sigmaUpdate] = computeOptimalUpdate(A, res, rho)
- rhoArho = scalarProduct(rho, A*rho);
- if max(abs(rho)) < eps
- lambda = 1;
- else
- lambda = scalarProduct(res, rho) ./ rhoArho;
- end
- sigmaUpdate = lambda.*rho;
- etaUpdate = 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
-
-function a = scalarProduct(x, y)
- a = sum(x.*y, 1);
-end
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
index fa30db22d2f6701c3fc1c43c78c4b31d43c63937..9bc2cf95231e39f832f7b21faad1791314702795 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/chooseIterativeSolver.m
@@ -55,7 +55,7 @@ switch class
if order == 1
solver = OptimalVcycleMultigridSolver(P1JacobiSmoother(fes, blf, P));
else
- solver = LocalMgLowOrderVcycle(fes, blf);
+ solver = OptimalVcycleMultigridSolver(JacobiLowOrderCascadeSmoother(fes, blf));
end
case "highOrderVcycle"
if order == 1
diff --git a/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m b/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m
new file mode 100644
index 0000000000000000000000000000000000000000..ad5dd1e038e06072b08b35a03d78f9514ef3af4b
--- /dev/null
+++ b/lib/solvers/smoothers/JacobiLowOrderCascadeSmoother.m
@@ -0,0 +1,88 @@
+% JacobiLowOrderCascadeSmoother (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 JacobiLowOrderCascadeSmoother < MultilevelSmoother
+ properties (Access=protected)
+ P
+ loFes
+ p1Smoother
+ hoSmoother
+ inclusionMatrix
+ intergridMatrix
+ freeVertices
+ freeVerticesOld
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = JacobiLowOrderCascadeSmoother(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
+ % intermediate level smoothers and transfer operators
+ obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
+ obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
+ 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
+
+ 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