From f3dd033e1bfad7a25a4c7a9e89aed53a84d48cc5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 9 Oct 2025 19:12:02 +0100 Subject: [PATCH 1/5] Quadrature elements constructed from non-Lagrange elements --- finat/quadrature.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/finat/quadrature.py b/finat/quadrature.py index 0272313c..11d3d2e1 100644 --- a/finat/quadrature.py +++ b/finat/quadrature.py @@ -26,6 +26,12 @@ def make_quadrature(ref_el, degree, scheme="default"): :arg degree: The degree of polynomial that the rule should integrate exactly. """ + + if hasattr(scheme, "dual_basis"): + Q, point_set = scheme.dual_basis + weights = numpy.ones((len(point_set.points),)) + return QuadratureRule(point_set, weights, ref_el=ref_el) + if ref_el.get_shape() == TENSORPRODUCT: try: degree = tuple(degree) From 4f2ed91ecb2337744eabdac01a5e0682edc5948d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 16 Oct 2025 10:27:18 +0100 Subject: [PATCH 2/5] Set weights to nan --- finat/quadrature.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/finat/quadrature.py b/finat/quadrature.py index 11d3d2e1..88b6c7ec 100644 --- a/finat/quadrature.py +++ b/finat/quadrature.py @@ -29,7 +29,8 @@ def make_quadrature(ref_el, degree, scheme="default"): if hasattr(scheme, "dual_basis"): Q, point_set = scheme.dual_basis - weights = numpy.ones((len(point_set.points),)) + # point_set does not define an integration rule + weights = numpy.full((len(point_set.points),), numpy.nan) return QuadratureRule(point_set, weights, ref_el=ref_el) if ref_el.get_shape() == TENSORPRODUCT: From df239521f33728df092623238c432cee27febe19 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 18 Dec 2025 09:25:28 +0000 Subject: [PATCH 3/5] Update finat/quadrature.py --- finat/quadrature.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/finat/quadrature.py b/finat/quadrature.py index 88b6c7ec..ed16d125 100644 --- a/finat/quadrature.py +++ b/finat/quadrature.py @@ -27,11 +27,9 @@ def make_quadrature(ref_el, degree, scheme="default"): integrate exactly. """ - if hasattr(scheme, "dual_basis"): - Q, point_set = scheme.dual_basis - # point_set does not define an integration rule - weights = numpy.full((len(point_set.points),), numpy.nan) - return QuadratureRule(point_set, weights, ref_el=ref_el) + if instance(scheme, QuadratureRule): + assert scheme.ref_el is ref_el + return scheme if ref_el.get_shape() == TENSORPRODUCT: try: From 9a33f5a9d47c0e872489d8b1f67c9a69f8168bde Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 18 Dec 2025 09:37:51 +0000 Subject: [PATCH 4/5] Fixup --- finat/quadrature.py | 5 ----- finat/quadrature_element.py | 6 +++++- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/finat/quadrature.py b/finat/quadrature.py index ed16d125..0272313c 100644 --- a/finat/quadrature.py +++ b/finat/quadrature.py @@ -26,11 +26,6 @@ def make_quadrature(ref_el, degree, scheme="default"): :arg degree: The degree of polynomial that the rule should integrate exactly. """ - - if instance(scheme, QuadratureRule): - assert scheme.ref_el is ref_el - return scheme - if ref_el.get_shape() == TENSORPRODUCT: try: degree = tuple(degree) diff --git a/finat/quadrature_element.py b/finat/quadrature_element.py index d56dd024..a6d0f308 100644 --- a/finat/quadrature_element.py +++ b/finat/quadrature_element.py @@ -31,7 +31,11 @@ def make_quadrature_element(fiat_ref_cell, degree, scheme="default", codim=0): else: rule_ref_cell = fiat_ref_cell - rule = make_quadrature(rule_ref_cell, degree, scheme=scheme) + if isinstance(scheme, AbstractQuadratureRule): + rule = scheme + else: + rule = make_quadrature(rule_ref_cell, degree, scheme=scheme) + return QuadratureElement(fiat_ref_cell, rule) From 9ee53b395cbb1060d6906aa03bb19bd03df62265 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 18 Dec 2025 12:46:48 +0000 Subject: [PATCH 5/5] test --- finat/quadrature_element.py | 1 + test/finat/test_quadrature_element.py | 23 +++++++++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 test/finat/test_quadrature_element.py diff --git a/finat/quadrature_element.py b/finat/quadrature_element.py index a6d0f308..818a70dd 100644 --- a/finat/quadrature_element.py +++ b/finat/quadrature_element.py @@ -33,6 +33,7 @@ def make_quadrature_element(fiat_ref_cell, degree, scheme="default", codim=0): if isinstance(scheme, AbstractQuadratureRule): rule = scheme + assert rule.ref_el == rule_ref_cell else: rule = make_quadrature(rule_ref_cell, degree, scheme=scheme) diff --git a/test/finat/test_quadrature_element.py b/test/finat/test_quadrature_element.py new file mode 100644 index 00000000..f6834b08 --- /dev/null +++ b/test/finat/test_quadrature_element.py @@ -0,0 +1,23 @@ +import pytest + +from FIAT import ufc_cell +from finat.quadrature import make_quadrature +from finat.quadrature_element import make_quadrature_element + + +@pytest.fixture(params=["interval", "triangle", "interval * interval", "triangle * interval"]) +def cell(request): + return ufc_cell(request.param) + + +def test_create_from_quadrature(cell): + degree = 4 + scheme = "default" + fe1 = make_quadrature_element(cell, degree, scheme=scheme) + + quadrature = make_quadrature(cell, degree, scheme=scheme) + fe2 = make_quadrature_element(cell, degree, scheme=quadrature) + + Q1, ps1 = fe1.dual_basis + Q2, ps2 = fe2.dual_basis + assert ps1.almost_equal(ps2)