diff --git a/lib/assembly/@BilinearForm/BilinearForm.m b/lib/assembly/@BilinearForm/BilinearForm.m
index 1221dd7053b48b1cea9af6cce316a872aa0add83..5387cef70ff48a924ebfca72f8ce7e23814520a9 100644
--- a/lib/assembly/@BilinearForm/BilinearForm.m
+++ b/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
diff --git a/lib/assembly/@BilinearForm/assemble.m b/lib/assembly/@BilinearForm/assemble.m
index 7a33f904e97f751d8983841258ef6356cad3cf64..cf56fa2b1c8abeefdadafcc345d9abe79c3b006e 100644
--- a/lib/assembly/@BilinearForm/assemble.m
+++ b/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
diff --git a/lib/assembly/@BilinearForm/computeVolumeData.m b/lib/assembly/@BilinearForm/computeVolumeData.m
index d13e55444301092bdfbaea550391ba11608ac4c5..3ff54fe51bae08ee2c498d8b9fbf05dce43c2200 100644
--- a/lib/assembly/@BilinearForm/computeVolumeData.m
+++ b/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
diff --git a/lib/assembly/@LinearForm/LinearForm.m b/lib/assembly/@LinearForm/LinearForm.m
index bbc44345e89e7594ec46e75a46baaad8824c4801..02e33ad5de8017998483fb554848153c9d7c9d89 100644
--- a/lib/assembly/@LinearForm/LinearForm.m
+++ b/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
diff --git a/lib/assembly/@LinearForm/assemble.m b/lib/assembly/@LinearForm/assemble.m
index 9c027b6caf887dd7867ec70cbce826fc44c29465..af9f28a0ce2b3ec69a36453318cd67dff2aef1d2 100644
--- a/lib/assembly/@LinearForm/assemble.m
+++ b/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
diff --git a/lib/quadrature/quadratureRules/QuadratureRule.m b/lib/quadrature/quadratureRules/QuadratureRule.m
index 38d2340340b1e091e0ce185a3ff0d45e65f630ad..ee69a439465fceca3d349bbd8baa3e8bf3bfe0ec 100644
--- a/lib/quadrature/quadratureRules/QuadratureRule.m
+++ b/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
 
diff --git a/lib/quadrature/quadratureRules/UnspecifiedQR.m b/lib/quadrature/quadratureRules/UnspecifiedQR.m
new file mode 100644
index 0000000000000000000000000000000000000000..9979bc84b7d7eb2a17ae0a6d0f1591b267c94aea
--- /dev/null
+++ b/lib/quadrature/quadratureRules/UnspecifiedQR.m
@@ -0,0 +1,14 @@
+% 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