diff --git a/lib/solvers/PcgSolver.m b/lib/solvers/PcgSolver.m
index 3a21ecc58737624493b4679cdb3dbfe9f211c153..1e8cd88c10f3074eaea6cf1c8eaad9682989a789 100644
--- a/lib/solvers/PcgSolver.m
+++ b/lib/solvers/PcgSolver.m
@@ -10,16 +10,30 @@ classdef PcgSolver < IterativeSolver
     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)
@@ -27,7 +41,7 @@ classdef PcgSolver < IterativeSolver
             
             % initialize residual & search direction
             obj.residual = b - obj.A*obj.x;
-            obj.Cresidual = obj.preconditionAction(obj.residual);
+            obj.Cresidual = obj.C.apply(obj.residual);
             obj.searchDirection = obj.Cresidual;
             obj.residualCNorm = sum(obj.residual.*obj.Cresidual, 1);
             obj.normb = sqrt(sum(b.^2, 1));
@@ -62,7 +76,7 @@ classdef PcgSolver < IterativeSolver
             
             % update residual
             obj.residual(:,idx) = obj.residual(:,idx) - alpha .* AsearchDirection;
-            obj.Cresidual(:,idx) = obj.preconditionAction(obj.residual(:,idx));
+            obj.Cresidual(:,idx) = obj.C.apply(obj.residual(:,idx));
             residualCNormOld = obj.residualCNorm(:,idx);
             obj.residualCNorm(:,idx) = sum(obj.residual(:,idx).*obj.Cresidual(:,idx), 1);
             
@@ -71,9 +85,4 @@ classdef PcgSolver < IterativeSolver
             obj.searchDirection(:,idx) = obj.Cresidual(:,idx) + beta .* obj.searchDirection(:,idx);
         end
     end
-    
-    %% abstract preconditioning
-    methods (Abstract, Access=public)
-        preconditionAction(obj, x)
-    end
 end
\ No newline at end of file
diff --git a/lib/solvers/PcgSolverTEMP.m b/lib/solvers/PcgSolverTEMP.m
deleted file mode 100644
index 53133bc01e43c65ae74528c0dfa167cc58355995..0000000000000000000000000000000000000000
--- a/lib/solvers/PcgSolverTEMP.m
+++ /dev/null
@@ -1,88 +0,0 @@
-% PcgSolver (abstract subclass of IterativeSolver) Solves linear equations
-%   iteratively by a preconditioned conjugate gradient (PCG) method.
-%
-% See also: IterativeSolver
-
-classdef PcgSolverTEMP < IterativeSolver
-    %% properties
-    properties (SetAccess=protected, GetAccess=public)
-        residualCNorm
-    end
-    
-    properties (Access=protected)
-        C
-        residual
-        Cresidual
-        searchDirection
-        normb
-    end
-
-    %% constructor
-    methods (Access=public)
-        function obj = PcgSolverTEMP(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 = sum(obj.residual.*obj.Cresidual, 1);
-            obj.normb = sqrt(sum(b.^2, 1));
-        end
-    end
-    
-    %% implement abstract superclass methods
-    methods (Access=public)
-        function tf = isConverged(obj)
-            tf = ((obj.iterationCount >= obj.maxIter) ...
-                        | (sqrt(obj.residualCNorm) < obj.tol) ...
-                        | (sqrt(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);
-            if sum(obj.searchDirection(:,idx).*AsearchDirection, 1) < eps
-                alpha = 1;
-            else
-                alpha = obj.residualCNorm(:,idx) ./ sum(obj.searchDirection(:,idx).*AsearchDirection, 1);
-            end
-            obj.x(:,idx) = obj.x(:,idx) + alpha .* obj.searchDirection(:,idx);
-
-            % DEBUG:
-            % disp(['alpha = ', num2str(alpha)])
-            
-            % 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) = sum(obj.residual(:,idx).*obj.Cresidual(:,idx), 1);
-            
-            % update search direction
-            beta = obj.residualCNorm(:,idx) ./ residualCNormOld;
-            obj.searchDirection(:,idx) = obj.Cresidual(:,idx) + beta .* obj.searchDirection(:,idx);
-        end
-    end
-end
\ No newline at end of file
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/chooseIterativeSolver.m
index 02e51cf7c5eafa6544617b44d5667a253fd79a2f..02d735c3b8fb17f82f0f0c2bae9c2a5658f8383a 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/chooseIterativeSolver.m
@@ -18,7 +18,7 @@ P = Prolongation.chooseFor(fes);
 switch class
     % non-preconditioned CG
     case "cg"
-        solver = PcgSolverTEMP(NoPreconditioner());
+        solver = PcgSolver(NoPreconditioner());
         
     % preconditioned CG family
     case "pcg"
@@ -46,7 +46,7 @@ switch class
             otherwise
                 error('No PCG variant %s!', variant)
         end
-        solver = PcgSolverTEMP(preconditioner);
+        solver = PcgSolver(preconditioner);
         
     % geometric multigrid family
     case "multigrid"