From 134e156b1ea6fbe85a6ab051cdb9d251c1f6743c Mon Sep 17 00:00:00 2001
From: Michael Innerberger <michael.innerberger@asc.tuwien.ac.at>
Date: Sat, 7 Oct 2023 05:41:21 -0400
Subject: [PATCH] Change way default values for quadrature rules are handled
---
lib/assembly/@BilinearForm/BilinearForm.m | 8 ++++----
lib/assembly/@BilinearForm/assemble.m | 3 ++-
lib/assembly/@BilinearForm/computeVolumeData.m | 10 +++++++---
lib/assembly/@LinearForm/LinearForm.m | 8 ++++----
lib/assembly/@LinearForm/assemble.m | 13 +++++++++----
.../quadratureRules/QuadratureRule.m | 18 ++++++++++++++++++
lib/quadrature/quadratureRules/UnspecifiedQR.m | 14 ++++++++++++++
7 files changed, 58 insertions(+), 16 deletions(-)
create mode 100644 lib/quadrature/quadratureRules/UnspecifiedQR.m
diff --git a/lib/assembly/@BilinearForm/BilinearForm.m b/lib/assembly/@BilinearForm/BilinearForm.m
index 1221dd7..5387cef 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 7a33f90..cf56fa2 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 d13e554..3ff54fe 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 bbc4434..02e33ad 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 9c027b6..af9f28a 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 38d2340..ee69a43 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 0000000..9979bc8
--- /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
--
GitLab