Change way default values for quadrature rules are handled

lib/assembly/@BilinearForm/BilinearForm.m

@@ -9,10 +9,10 @@ classdef BilinearForm < handle
         b {mustBeEvaluableOrEmpty} = []     % Convection vector field
         c {mustBeEvaluableOrEmpty} = []     % Reaction coefficient
         robin {mustBeEvaluableOrEmpty} = [] % Robin coefficient
-        qra (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
-        qrb (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
-        qrc (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
-        qrRobin (1,1) QuadratureRule = QuadratureRule.ofOrder(1, '1D')
+        qra (1,1) QuadratureRule = UnspecifiedQR
+        qrb (1,1) QuadratureRule = UnspecifiedQR
+        qrc (1,1) QuadratureRule = UnspecifiedQR
+        qrRobin (1,1) QuadratureRule = UnspecifiedQR
     end
     
     %% methods

lib/assembly/@BilinearForm/assemble.m

@@ -25,12 +25,13 @@ end
 %% integrate and add boundary data
 % Robin data
 if ~isempty(obj.robin)
+    qrRobin = QuadratureRule.ifEmptyOfOrder(obj.qrRobin, 2*fes.finiteElement.order, "1D");
     phi = TestFunction(fes);
     f = CompositeFunction(@BilinearForm.robinPart, obj.robin, phi);
     idx = getCombinedBndEdges(fes.mesh, fes.bnd.robin);
     [I, J] = getLocalDofs(size(dofs.edge2Dofs, Dim.Vector));
     mat = mat + sparse(dofs.edge2Dofs(I,idx), dofs.edge2Dofs(J,idx), ...
-        integrateEdge(f, obj.qrRobin, idx), dofs.nDofs, dofs.nDofs);
+        integrateEdge(f, qrRobin, idx), dofs.nDofs, dofs.nDofs);
 end
 
 end

lib/assembly/@BilinearForm/computeVolumeData.m

@@ -19,23 +19,27 @@ end
 data = 0;
 phi = TestFunction(fes);
 Dphi = TestFunctionGradient(fes);
+p = fes.finiteElement.order;
 
 % diffusion
 if ~isempty(obj.a)
+    qr = QuadratureRule.ifEmptyOfOrder(obj.qra, 2*(p-1));
     f = CompositeFunction(@BilinearForm.diffusionPart, obj.a, Dphi);
-    data = data + integrateElement(f, obj.qra);
+    data = data + integrateElement(f, qr);
 end
 
 % convection
 if ~isempty(obj.b)
+    qr = QuadratureRule.ifEmptyOfOrder(obj.qrb, 2*p-1);
     f = CompositeFunction(@BilinearForm.convectionPart, obj.b, Dphi, phi);
-    data = data + integrateElement(f, obj.qrb);
+    data = data + integrateElement(f, qr);
 end
 
 % reaction
 if ~isempty(obj.c)
+    qr = QuadratureRule.ifEmptyOfOrder(obj.qrc, 2*p);
     f = CompositeFunction(@BilinearForm.reactionPart, obj.c, phi);
-    data = data + integrateElement(f, obj.qrc);
+    data = data + integrateElement(f, qr);
 end
 
 end

lib/assembly/@LinearForm/LinearForm.m

@@ -9,10 +9,10 @@ classdef LinearForm < handle
         fvec {mustBeEvaluableOrEmpty} = []      % Generalized volume load
         robin {mustBeEvaluableOrEmpty} = []     % Robin boundary load
         neumann {mustBeEvaluableOrEmpty} = []   % Neumann boundary load
-        qrf (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
-        qrfvec (1,1) QuadratureRule = QuadratureRule.ofOrder(1)
-        qrRobin (1,1) QuadratureRule = QuadratureRule.ofOrder(1, '1D')
-        qrNeumann (1,1) QuadratureRule = QuadratureRule.ofOrder(1, '1D')
+        qrf (1,1) QuadratureRule = UnspecifiedQR
+        qrfvec (1,1) QuadratureRule = UnspecifiedQR
+        qrRobin (1,1) QuadratureRule = UnspecifiedQR
+        qrNeumann (1,1) QuadratureRule = UnspecifiedQR
     end
     
     %% methods

lib/assembly/@LinearForm/assemble.m

@@ -19,17 +19,20 @@ end
 vec = 0;
 phi = TestFunction(fes);
 Dphi = TestFunctionGradient(fes);
+p = fes.finiteElement.order;
 
 % scalar rhs
 if ~isempty(obj.f)
+    qr = QuadratureRule.ifEmptyOfOrder(obj.qrf, p);
     f = CompositeFunction(@LinearForm.scalarPart, obj.f, phi);
-    vec = vec + integrateElement(f, obj.qrf);
+    vec = vec + integrateElement(f, qr);
 end
 
 % vector rhs
 if ~isempty(obj.fvec)
+    qr = QuadratureRule.ifEmptyOfOrder(obj.qrfvec, p-1);
     f = CompositeFunction(@LinearForm.vectorPart, obj.fvec, Dphi);
-    vec = vec + integrateElement(f, obj.qrfvec);
+    vec = vec + integrateElement(f, qr);
 end
 
 %% construct rhs vector from element data
@@ -39,18 +42,20 @@ vec = accumarray(dofs.element2Dofs(:), vec(:), [dofs.nDofs, 1]);
 %% integrate and add boundary data
 % Robin data
 if ~isempty(obj.robin)
+    qr = QuadratureRule.ifEmptyOfOrder(obj.qrRobin, p, "1D");
     f = CompositeFunction(@LinearForm.robinPart, obj.robin, phi);
     idx = getCombinedBndEdges(fes.mesh, fes.bnd.robin);
-    edgeData = integrateEdge(f, obj.qrRobin, idx);
+    edgeData = integrateEdge(f, qr, idx);
     edge2Dofs = dofs.edge2Dofs(:,idx);
     vec = vec + accumarray(edge2Dofs(:), edgeData(:), [dofs.nDofs, 1]);
 end
 
 % Neumann data
 if ~isempty(obj.neumann)
+    qr = QuadratureRule.ifEmptyOfOrder(obj.qrNeumann, p-1, "1D");
     f = CompositeFunction(@LinearForm.neumannPart, obj.neumann, phi);
     idx = getCombinedBndEdges(fes.mesh, fes.bnd.neumann);
-    edgeData = integrateEdge(f, obj.qrNeumann, idx);
+    edgeData = integrateEdge(f, qr, idx);
     edge2Dofs = dofs.edge2Dofs(:,idx);
     vec = vec + accumarray(edge2Dofs(:), edgeData(:), [dofs.nDofs, 1]);
 end

lib/quadrature/quadratureRules/QuadratureRule.m

@@ -39,6 +39,7 @@ classdef QuadratureRule
     end
     
     methods(Static)
+        % Return quadrature rule of specified order and dimension
         function obj = ofOrder(order, dim)
             arguments
                 order {mustBePositive(order)}
@@ -56,6 +57,23 @@ classdef QuadratureRule
             
             obj = QuadratureRule(bary, weights);
         end
+
+        % Return given quadrature rule or create new one if given rule is
+        % an UnspecifiedQR
+        function qr = ifEmptyOfOrder(qrOrUnspecified, order, dim)
+            arguments
+                qrOrUnspecified (1,1) QuadratureRule
+                order (1,1) double
+                dim {mustBeTextScalar, mustBeMember(dim, {'1D', '2D'})} = '2D'
+            end
+
+            if isa(qrOrUnspecified, 'UnspecifiedQR')
+                order = max(order, 1);
+                qr = QuadratureRule.ofOrder(order, dim);
+            else
+                qr = qrOrUnspecified;
+            end
+        end
     end
 end
 

lib/quadrature/quadratureRules/UnspecifiedQR.m

+% UnspecifiedQR Null object as placeholder for not-set quadrature rule.
+
+classdef UnspecifiedQR < QuadratureRule
+    %% methods
+    methods (Access=public)
+        function obj = UnspecifiedQR()
+            obj = obj@QuadratureRule(Barycentric1D([0; 1]), 1);
+            obj.bary = [];
+            obj.weights = [];
+            obj.nNodes = 0;
+            obj.dim = NaN;
+        end
+    end
+end
\ No newline at end of file