diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
index 877bda21be9437ae223500a29d24fa4b25b7aec9..f353744f765d7017efa7c02823f5d3f8c96a3b8a 100644
--- a/examples/algebraicSolver.m
+++ b/examples/algebraicSolver.m
@@ -24,10 +24,7 @@ for p = 1:pmax
     lf = LinearForm();
 
     blf.a = Constant(mesh, 1);
-    blf.qra = QuadratureRule.ofOrder(max(2*p-2, 1));
-
     lf.f = Constant(mesh, 1);
-    lf.qrf = QuadratureRule.ofOrder(2*p);
 
     %% set up solver and operator for nested iteration
     % choose the iterative solver of choice: multigrid (with the variants
diff --git a/examples/convergenceRates.m b/examples/convergenceRates.m
index f8b6455f400b32a5dcf74fa926cfa62d4284fed5..d054129fbbc44f028d3f3ff90f5f06e06f5e69e8 100644
--- a/examples/convergenceRates.m
+++ b/examples/convergenceRates.m
@@ -22,9 +22,9 @@ for p = 1:pmax
     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');
 
diff --git a/examples/goafemIterativeSolver.m b/examples/goafemIterativeSolver.m
index 357b256be0535e36cdfd97636a33d70296c106e2..aed69d04e8ac2983c895c06b27299ea46a611128 100644
--- a/examples/goafemIterativeSolver.m
+++ b/examples/goafemIterativeSolver.m
@@ -26,7 +26,6 @@ for p = 1:pmax
     % dual:   -\Delta z = "lorentzian peak 2"
     blf = BilinearForm();
     blf.a = Constant(mesh, 1);
-    blf.qra = QuadratureRule.ofOrder(max(2*p-2, 1));
 
     lfF = LinearForm();
     lfF.f = MeshFunction(mesh, @(x) lorentzian(x, [0.7;0.7], 1e-1));
diff --git a/examples/goalOrientedAFEM.m b/examples/goalOrientedAFEM.m
index 8cc4bc6b5fa244ee2aa2779e346adcd9beddcfa2..a07806216b2b6a115b00a751585e3797c48ba33d 100644
--- a/examples/goalOrientedAFEM.m
+++ b/examples/goalOrientedAFEM.m
@@ -22,21 +22,19 @@ for p = 1:pmax
     z = FeFunction(fes);
     
     %% set problem data given in [Mommer, Stevenson; 2009]
+	% 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();
     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();
     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 = LoMeshProlongation(ncFes);
diff --git a/examples/iterativeLinearization.m b/examples/iterativeLinearization.m
index e4bf7d3b7c2dfe85e1f7485ab22b24c5f5aa97c6..57a6f1782c0382497796e65eda21b310ba3677e7 100644
--- a/examples/iterativeLinearization.m
+++ b/examples/iterativeLinearization.m
@@ -116,8 +116,6 @@ function [blf, lf] = setProblemData(mesh, 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)
diff --git a/experiments/p_laplace.m b/experiments/p_laplace.m
index 9bee4bdbb032684aa188423747fe374f545fd6f1..16e81435101ff679cd305a62701dcc28128c05f7 100644
--- a/experiments/p_laplace.m
+++ b/experiments/p_laplace.m
@@ -113,8 +113,6 @@ function [blf, lf] = setProblemData(mesh, fes, Du, g, linearization, p)
             lf.f = g;
             lf.fvec = CompositeFunction(@(t) -mu(vectorProduct(t, t),p) .* t, Du);
     end
-
-    [blf.qra, lf.qrf, lf.qrfvec] = deal(QuadratureRule.ofOrder(1));
 end
 
 function updateDataU(u, deltaZ, A, F, linearization)
diff --git a/runAdditiveSchwarzPreconditioner.m b/runAdditiveSchwarzPreconditioner.m
index 971ffd8d1f8979c73c28c1216fafa76b00187343..a0dbc376bbeba500ba22ff0476058fae1069bdce 100644
--- a/runAdditiveSchwarzPreconditioner.m
+++ b/runAdditiveSchwarzPreconditioner.m
@@ -28,7 +28,6 @@ function leveldata = runAdditiveSchwarzPreconditioner(pMax)
         blf = BilinearForm();
         lf = LinearForm();
         blf.a = Constant(mesh, 1);
-        blf.qra = QuadratureRule.ofOrder(max(2*(p-1), 1));
         lf.f = Constant(mesh, 1);
     
         % Choose iterative solver
diff --git a/runIterativeSolver.m b/runIterativeSolver.m
index ce2bcf8b6d8cd180066937d7c5843e07e05907f6..a14ceae87e6bf7a7b550c59b7b2b0b33b338c0ec 100644
--- a/runIterativeSolver.m
+++ b/runIterativeSolver.m
@@ -24,7 +24,6 @@ function leveldata = runIterativeSolver(maxNiter)
     blf = BilinearForm();
     lf = LinearForm();
     blf.a = Constant(mesh, 1);
-    blf.qra = QuadratureRule.ofOrder(max(2*(p-1), 1));
     lf.f = Constant(mesh, 1);
 
     % Choose iterative solver
diff --git a/tests/TestIterativeSolver.m b/tests/TestIterativeSolver.m
index 6a32e6af2f72e24d59a80be2ceb5d94690187682..11964911c3bc50586be4a40a1b716d06cf1dc175 100644
--- a/tests/TestIterativeSolver.m
+++ b/tests/TestIterativeSolver.m
@@ -74,9 +74,7 @@ methods (Access=private)
         blf = BilinearForm();
         lf = LinearForm();
         blf.a = Constant(mesh, 1);
-        blf.qra = QuadratureRule.ofOrder(max(2*(p-1), 1));
         blf.c = Constant(mesh, 1);
-        blf.qrc = QuadratureRule.ofOrder(2*p);
         lf.f = Constant(mesh, 1);
     end
     
diff --git a/tests/TestPatchwiseMatrix.m b/tests/TestPatchwiseMatrix.m
index d44b7b64dd3a32b36406483f4e93764cc8bdd63b..206b6a6351075c53ed9b1797c5e783ed868cbaac 100644
--- a/tests/TestPatchwiseMatrix.m
+++ b/tests/TestPatchwiseMatrix.m
@@ -49,7 +49,6 @@ methods (Test)
         fes = FeSpace(mesh, HigherOrderH1Fe(p), 'dirichlet', [], 'neumann', ':');
         blf = BilinearForm();
         blf.a = Constant(mesh, 1);
-        blf.qra = QuadratureRule.ofOrder(max(2*(p-1),1));
         blf.c = MeshFunction(mesh, @(x) 1+x(1,:,:));
         blf.qrc = QuadratureRule.ofOrder(2*p+1);
         A = assemble(blf, fes);