From e2b6459b48696d85a81a25bc01e965f9329475f1 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <michael.innerberger@asc.tuwien.ac.at>
Date: Sat, 7 Oct 2023 06:06:47 -0400
Subject: [PATCH] Remove explicit quadrature rule where applicable

---
 examples/algebraicSolver.m         | 3 ---
 examples/convergenceRates.m        | 4 ++--
 examples/goafemIterativeSolver.m   | 1 -
 examples/goalOrientedAFEM.m        | 4 +---
 examples/iterativeLinearization.m  | 2 --
 experiments/p_laplace.m            | 2 --
 runAdditiveSchwarzPreconditioner.m | 1 -
 runIterativeSolver.m               | 1 -
 tests/TestIterativeSolver.m        | 2 --
 tests/TestPatchwiseMatrix.m        | 1 -
 10 files changed, 3 insertions(+), 18 deletions(-)

diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
index 877bda2..f353744 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 f8b6455..d054129 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 357b256..aed69d0 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 8cc4bc6..a078062 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 e4bf7d3..57a6f17 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 9bee4bd..16e8143 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 971ffd8..a0dbc37 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 ce2bcf8..a14ceae 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 6a32e6a..1196491 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 d44b7b6..206b6a6 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);
-- 
GitLab