From 657321796eeda75085e49048cc759f70d587cc5b Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 19 Aug 2024 13:12:40 +0100 Subject: [PATCH 01/35] First draft of expanding set size to match when unioning --- FIAT/polynomial_set.py | 38 +++++++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 732d03892..a3232eba1 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -169,7 +169,9 @@ def polynomial_set_union_normalized(A, B): not contain any of the same members of the set, as we construct a span via SVD. """ - new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) + # new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) + assert A.get_reference_element() == B.get_reference_element() + new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B) func_shape = new_coeffs.shape[1:] if len(func_shape) == 1: (u, sig, vt) = numpy.linalg.svd(new_coeffs) @@ -192,6 +194,40 @@ def polynomial_set_union_normalized(A, B): coeffs) +def construct_new_coeffs(ref_el, A, B): + # Constructs new coefficients for the set union of A and B + # If A and B do not have the same degree the smaller one + # is extended to match the larger + + sd = ref_el.get_spatial_dimension() + if A.degree != B.degree: + higher = A if A.degree > B.degree else B + lower = B if A.degree > B.degree else A + + dimHigher = expansions.polynomial_dimension(ref_el, higher.degree) + dimLower = expansions.polynomial_dimension(ref_el, lower.degree) + + if (dimLower == len(list(lower.coeffs))): + lower_indices = list(chain(*(range(i * dimHigher, i * dimHigher + dimLower) for i in range(sd)))) + embedded = higher.take(lower_indices) + embedded_coeffs = embedded.coeffs + else: + # if lower space not complete take a different approach + embedded_coeffs = [] + diff = dimHigher - dimLower + for coeff in lower.coeffs: + new_coeff = [] + for row in coeff: + new_coeff.append(numpy.append(row, [0 for i in range(diff)])) + embedded_coeffs.append(new_coeff) + + new_coeffs = numpy.array(list(embedded_coeffs) + list(higher.coeffs)) + else: + new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) + + return new_coeffs + + class ONSymTensorPolynomialSet(PolynomialSet): """Constructs an orthonormal basis for symmetric-tensor-valued polynomials on a reference element. From dc5188a822e684c9466111645dcf7ca1fd9d1dc7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 19 Aug 2024 14:06:49 +0100 Subject: [PATCH 02/35] Change condition of extension --- FIAT/polynomial_set.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index a3232eba1..0b64c3cf7 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -200,7 +200,7 @@ def construct_new_coeffs(ref_el, A, B): # is extended to match the larger sd = ref_el.get_spatial_dimension() - if A.degree != B.degree: + if A.get_embedded_degree() != B.get_embedded_degree(): higher = A if A.degree > B.degree else B lower = B if A.degree > B.degree else A From c75efef8dc904005213ece3877a9f04835d3ec0e Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 20 Aug 2024 14:53:13 +0100 Subject: [PATCH 03/35] Adapt new coeff construction for vec sets --- FIAT/polynomial_set.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 0b64c3cf7..27292c171 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -204,23 +204,28 @@ def construct_new_coeffs(ref_el, A, B): higher = A if A.degree > B.degree else B lower = B if A.degree > B.degree else A - dimHigher = expansions.polynomial_dimension(ref_el, higher.degree) - dimLower = expansions.polynomial_dimension(ref_el, lower.degree) - - if (dimLower == len(list(lower.coeffs))): - lower_indices = list(chain(*(range(i * dimHigher, i * dimHigher + dimLower) for i in range(sd)))) - embedded = higher.take(lower_indices) - embedded_coeffs = embedded.coeffs - else: - # if lower space not complete take a different approach - embedded_coeffs = [] - diff = dimHigher - dimLower - for coeff in lower.coeffs: + # dimHigher = expansions.polynomial_dimension(ref_el, higher.degree) + # dimLower = expansions.polynomial_dimension(ref_el, lower.degree) + + try: + sd = lower.get_shape()[0] + except IndexError: + sd = 1 + embedded_coeffs = [] + diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1] + for coeff in lower.coeffs: + if sd > 1: new_coeff = [] for row in coeff: new_coeff.append(numpy.append(row, [0 for i in range(diff)])) embedded_coeffs.append(new_coeff) + else: + embedded_coeffs.append(numpy.append(coeff, [0 for i in range(diff)])) + embedded_coeffs = numpy.array(embedded_coeffs) + # print("embedded", embedded_coeffs.shape) + # print("higher", higher.coeffs.shape) + # print("lower", lower.coeffs.shape) new_coeffs = numpy.array(list(embedded_coeffs) + list(higher.coeffs)) else: new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) From dfd13475368d1ca3d935494747bbda6503f78961 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 21 Aug 2024 12:02:48 +0100 Subject: [PATCH 04/35] select larger of degrees --- FIAT/polynomial_set.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 27292c171..7553007dc 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -187,9 +187,11 @@ def polynomial_set_union_normalized(A, B): coeffs = numpy.reshape(vt[:num_sv], (num_sv,) + func_shape) + deg = max(A.get_degree(), B.get_degree()) + em_deg = max(A.get_embedded_degree(), B.get_embedded_degree()) return PolynomialSet(A.get_reference_element(), - A.get_degree(), - A.get_embedded_degree(), + deg, + em_deg, A.get_expansion_set(), coeffs) @@ -223,9 +225,6 @@ def construct_new_coeffs(ref_el, A, B): embedded_coeffs.append(numpy.append(coeff, [0 for i in range(diff)])) embedded_coeffs = numpy.array(embedded_coeffs) - # print("embedded", embedded_coeffs.shape) - # print("higher", higher.coeffs.shape) - # print("lower", lower.coeffs.shape) new_coeffs = numpy.array(list(embedded_coeffs) + list(higher.coeffs)) else: new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) From 21f1bf163a0ac185075d79da19ec4729c3a8a2b8 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 28 Aug 2024 14:35:07 +0100 Subject: [PATCH 05/35] Add orthonormal requirement --- FIAT/expansions.py | 2 +- FIAT/nedelec.py | 8 ++++---- FIAT/polynomial_set.py | 10 ++++------ FIAT/raviart_thomas.py | 6 +++--- 4 files changed, 12 insertions(+), 14 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index bf38f9795..bf2015b18 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -327,7 +327,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): sd = self.ref_el.get_spatial_dimension() # Always return 1 for n=0 to make regression tests pass - scale = 1.0 if n == 0 and len(self.affine_mappings) == 1 else self.get_scale(cell=cell) + scale = 1.0 if n == 0 and (self.scale is None) else self.get_scale(cell=cell) phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv, scale, variant=self.variant) if self.continuity == "C0": diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c05133022..ca135074c 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -22,7 +22,7 @@ def NedelecSpace2D(ref_el, degree): raise Exception("NedelecSpace2D requires 2d reference element") k = degree - 1 - vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,)) + vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal") dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1) dimPk = expansions.polynomial_dimension(ref_el, k) @@ -32,7 +32,7 @@ def NedelecSpace2D(ref_el, degree): for i in range(sd)))) vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) - Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1) + Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal") PkH = Pkp1.take(list(range(dimPkm1, dimPk))) Q = create_quadrature(ref_el, 2 * (k + 1)) @@ -43,8 +43,8 @@ def NedelecSpace2D(ref_el, degree): CrossX = numpy.dot(numpy.array([[0.0, 1.0], [-1.0, 0.0]]), Qpts.T) PkHCrossX_at_Qpts = PkH_at_Qpts[:, None, :] * CrossX[None, :, :] - PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), Pkp1_at_Qpts.T) - + PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), + Pkp1_at_Qpts.T) PkHcrossX = polynomial_set.PolynomialSet(ref_el, k + 1, k + 1, diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 7553007dc..4fb7720ac 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -169,10 +169,10 @@ def polynomial_set_union_normalized(A, B): not contain any of the same members of the set, as we construct a span via SVD. """ - # new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) assert A.get_reference_element() == B.get_reference_element() new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B) func_shape = new_coeffs.shape[1:] + if len(func_shape) == 1: (u, sig, vt) = numpy.linalg.svd(new_coeffs) num_sv = len([s for s in sig if abs(s) > 1.e-10]) @@ -182,7 +182,9 @@ def polynomial_set_union_normalized(A, B): new_shape1 = numpy.prod(func_shape) newshape = (new_shape0, new_shape1) nc = numpy.reshape(new_coeffs, newshape) + (u, sig, vt) = numpy.linalg.svd(nc, 1) + num_sv = len([s for s in sig if abs(s) > 1.e-10]) coeffs = numpy.reshape(vt[:num_sv], (num_sv,) + func_shape) @@ -206,13 +208,11 @@ def construct_new_coeffs(ref_el, A, B): higher = A if A.degree > B.degree else B lower = B if A.degree > B.degree else A - # dimHigher = expansions.polynomial_dimension(ref_el, higher.degree) - # dimLower = expansions.polynomial_dimension(ref_el, lower.degree) - try: sd = lower.get_shape()[0] except IndexError: sd = 1 + embedded_coeffs = [] diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1] for coeff in lower.coeffs: @@ -224,11 +224,9 @@ def construct_new_coeffs(ref_el, A, B): else: embedded_coeffs.append(numpy.append(coeff, [0 for i in range(diff)])) embedded_coeffs = numpy.array(embedded_coeffs) - new_coeffs = numpy.array(list(embedded_coeffs) + list(higher.coeffs)) else: new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) - return new_coeffs diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index f33ce937f..13de2e46f 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -20,7 +20,7 @@ def RTSpace(ref_el, degree): sd = ref_el.get_spatial_dimension() k = degree - 1 - vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,)) + vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal") dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1) dimPk = expansions.polynomial_dimension(ref_el, k) @@ -30,7 +30,7 @@ def RTSpace(ref_el, degree): for i in range(sd)))) vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) - Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1) + Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal") PkH = Pkp1.take(list(range(dimPkm1, dimPk))) Q = create_quadrature(ref_el, 2 * (k + 1)) @@ -39,12 +39,12 @@ def RTSpace(ref_el, degree): # have to work on this through "tabulate" interface # first, tabulate PkH at quadrature points PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] x = Qpts.T PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) - PkHx = polynomial_set.PolynomialSet(ref_el, k, k + 1, From 2b45bd575f08eb91b8b0a7d230a8d839033a8d54 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 9 Dec 2024 11:05:25 +0000 Subject: [PATCH 06/35] Move finat commits from old branch --- finat/__init__.py | 1 + finat/fiat_elements.py | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/finat/__init__.py b/finat/__init__.py index cad018212..4b71639a4 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -12,6 +12,7 @@ from .fiat_elements import GopalakrishnanLedererSchoberlSecondKind # noqa: F401 from .fiat_elements import FacetBubble # noqa: F401 from .fiat_elements import KongMulderVeldhuizen # noqa: F401 +from .fiat_elements import IndiaDefElement #noqa: F401 from .argyris import Argyris # noqa: F401 from .aw import ArnoldWinther # noqa: F401 diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 1f2081894..5876a2da9 100644 --- a/finat/fiat_elements.py +++ b/finat/fiat_elements.py @@ -480,3 +480,8 @@ def __init__(self, cell, degree, variant=None): class NedelecSecondKind(VectorFiatElement): def __init__(self, cell, degree, variant=None): super().__init__(FIAT.NedelecSecondKind(cell, degree, variant=variant)) + + +class IndiaDefElement(FiatElement): + def __init__(self, triple): + super(IndiaDefElement, self).__init__(triple.to_fiat_elem()) From 9b0568841ff51606bb52c9041f686bbc24a60ebf Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 9 Dec 2024 11:06:11 +0000 Subject: [PATCH 07/35] Lint --- finat/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/finat/__init__.py b/finat/__init__.py index 4b71639a4..7e70161a4 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -12,7 +12,7 @@ from .fiat_elements import GopalakrishnanLedererSchoberlSecondKind # noqa: F401 from .fiat_elements import FacetBubble # noqa: F401 from .fiat_elements import KongMulderVeldhuizen # noqa: F401 -from .fiat_elements import IndiaDefElement #noqa: F401 +from .fiat_elements import IndiaDefElement # noqa: F401 from .argyris import Argyris # noqa: F401 from .aw import ArnoldWinther # noqa: F401 From 8d1a0747402b9734807bb35420630e59f4f0aa63 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 6 Jan 2025 13:31:27 +0000 Subject: [PATCH 08/35] Refactor to deal with single source of truth --- finat/element_factory.py | 8 ++++++++ finat/ufl/__init__.py | 1 + finat/ufl/fuseelement.py | 42 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 51 insertions(+) create mode 100644 finat/ufl/fuseelement.py diff --git a/finat/element_factory.py b/finat/element_factory.py index 48db428d8..eead428ee 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -24,6 +24,7 @@ import finat import finat.ufl +from redefining_fe.cells import CellComplexToUFL as FuseCell import ufl from FIAT import ufc_cell @@ -112,6 +113,8 @@ def as_fiat_cell(cell): :arg cell: the :class:`ufl.Cell` to convert.""" if not isinstance(cell, ufl.AbstractCell): raise ValueError("Expecting a UFL Cell") + if isinstance(cell, FuseCell): + return cell.to_fiat() return ufc_cell(cell) @@ -302,6 +305,11 @@ def convert_restrictedelement(element, **kwargs): return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps +@convert.register(finat.ufl.FuseElement) +def convert_india_def(element, **kwargs): + return finat.fiat_elements.IndiaDefElement(element.triple), set() + + hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval) quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval) _cache = weakref.WeakKeyDictionary() diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py index 21a7d13d1..33373eb59 100644 --- a/finat/ufl/__init__.py +++ b/finat/ufl/__init__.py @@ -19,3 +19,4 @@ from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement # noqa: F401 from finat.ufl.restrictedelement import RestrictedElement # noqa: F401 from finat.ufl.tensorproductelement import TensorProductElement # noqa: F401 +from finat.ufl.fuseelement import FuseElement # noqa: F401 \ No newline at end of file diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py new file mode 100644 index 000000000..8b75fc676 --- /dev/null +++ b/finat/ufl/fuseelement.py @@ -0,0 +1,42 @@ +"""Element.""" +# -*- coding: utf-8 -*- +# Copyright (C) 2025 India Marsden +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +from finat.ufl.finiteelementbase import FiniteElementBase + +class FuseElement(FiniteElementBase): + """ + TODO: Need to deal with cases where value shape and reference value shape are different + """ + + def __init__(self, triple, cell=None): + self.triple = triple + if not cell: + cell = self.triple.cell.to_ufl() + + # this isn't really correct + degree = self.triple.spaces[0].degree() + super(FuseElement, self).__init__("IT", cell, degree, None, triple.get_value_shape()) + + def __repr__(self): + return "FiniteElement(%s, %s, (%s, %s, %s), %s)" % ( + repr(self.triple.DOFGenerator), repr(self.triple.cell), repr(self.triple.spaces[0]), repr(self.triple.spaces[1]), repr(self.triple.spaces[2]), "X") + + def __str__(self): + return "" % (self.triple.spaces[0], self.triple.cell) + + def mapping(self): + if str(self.sobolev_space) == "HCurl": + return "covariant Piola" + elif str(self.sobolev_space) == "HDiv": + return "contravariant Piola" + else: + return "identity" + + def sobolev_space(self): + return self.triple.spaces[1] + + def reconstruct(self, family=None, cell=None, degree=None, quad_scheme=None, variant=None): + return FuseElement(self.triple, cell=cell) From b72cd61da8b735ed14af0b2ed96d59c2c621266b Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 6 Jan 2025 15:49:12 +0000 Subject: [PATCH 09/35] Refactoring --- finat/fiat_elements.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 04e25089b..4ac4aa7ba 100644 --- a/finat/fiat_elements.py +++ b/finat/fiat_elements.py @@ -484,4 +484,4 @@ def __init__(self, cell, degree, variant=None): class IndiaDefElement(FiatElement): def __init__(self, triple): - super(IndiaDefElement, self).__init__(triple.to_fiat_elem()) + super(IndiaDefElement, self).__init__(triple.to_fiat()) From a89d057d51f6ea382182a1e2e1676b684c3803a7 Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Fri, 10 Jan 2025 11:17:20 +0000 Subject: [PATCH 10/35] Renaming fuse project (#125) Renaming + adding to install --- finat/element_factory.py | 2 +- finat/ufl/__init__.py | 2 +- finat/ufl/fuseelement.py | 3 +++ pyproject.toml | 1 + 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/finat/element_factory.py b/finat/element_factory.py index eead428ee..5adcca5f9 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -24,7 +24,7 @@ import finat import finat.ufl -from redefining_fe.cells import CellComplexToUFL as FuseCell +from fuse.cells import CellComplexToUFL as FuseCell import ufl from FIAT import ufc_cell diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py index 33373eb59..fd656db6e 100644 --- a/finat/ufl/__init__.py +++ b/finat/ufl/__init__.py @@ -19,4 +19,4 @@ from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement # noqa: F401 from finat.ufl.restrictedelement import RestrictedElement # noqa: F401 from finat.ufl.tensorproductelement import TensorProductElement # noqa: F401 -from finat.ufl.fuseelement import FuseElement # noqa: F401 \ No newline at end of file +from finat.ufl.fuseelement import FuseElement # noqa: F401 diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py index 8b75fc676..b83b26ba0 100644 --- a/finat/ufl/fuseelement.py +++ b/finat/ufl/fuseelement.py @@ -6,8 +6,11 @@ from finat.ufl.finiteelementbase import FiniteElementBase + class FuseElement(FiniteElementBase): """ + A finite element defined using FUSE. + TODO: Need to deal with cases where value shape and reference value shape are different """ diff --git a/pyproject.toml b/pyproject.toml index c224c1fbe..2d43b6dc5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,6 +12,7 @@ dependencies = [ "symengine", "sympy", "fenics-ufl @ git+https://github.com/firedrakeproject/ufl.git", + "fuse-element @ git+https://github.com/indiamai/fuse.git", ] requires-python = ">=3.10" authors = [ From ac81474f539c09adaca0b520325c68e37cf6f3fa Mon Sep 17 00:00:00 2001 From: India Marsden Date: Sat, 11 Jan 2025 11:54:04 +0000 Subject: [PATCH 11/35] Address Pablo's comments --- FIAT/polynomial_set.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index ce90dca14..7a1df1013 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -190,12 +190,18 @@ def polynomial_set_union_normalized(A, B): def construct_new_coeffs(ref_el, A, B): - # Constructs new coefficients for the set union of A and B - # If A and B do not have the same degree the smaller one - # is extended to match the larger + """ + Constructs new coefficients for the set union of A and B + If A and B are discontinuous and do not have the same degree the smaller one + is extended to match the larger. + + This does not handle the case that A and B have continuity but not the same degree. + """ sd = ref_el.get_spatial_dimension() - if A.get_embedded_degree() != B.get_embedded_degree(): + if A.get_expansion_set().continuity != B.get_expansion_set().continuity: + raise ValueError("Continuity of expansion sets does not match.") + if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None: higher = A if A.degree > B.degree else B lower = B if A.degree > B.degree else A @@ -215,9 +221,9 @@ def construct_new_coeffs(ref_el, A, B): else: embedded_coeffs.append(numpy.append(coeff, [0 for i in range(diff)])) embedded_coeffs = numpy.array(embedded_coeffs) - new_coeffs = numpy.array(list(embedded_coeffs) + list(higher.coeffs)) + new_coeffs = numpy.concatenate((embedded_coeffs, higher.coeffs), axis=0) else: - new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) + new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) return new_coeffs From 9420a6dbe6e94a10ed89aec50ffef6e9077e144b Mon Sep 17 00:00:00 2001 From: India Marsden Date: Sat, 11 Jan 2025 11:55:20 +0000 Subject: [PATCH 12/35] Lint --- FIAT/polynomial_set.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 7a1df1013..4f61968bd 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -194,7 +194,7 @@ def construct_new_coeffs(ref_el, A, B): Constructs new coefficients for the set union of A and B If A and B are discontinuous and do not have the same degree the smaller one is extended to match the larger. - + This does not handle the case that A and B have continuity but not the same degree. """ From be81de3f9297ac0b45390ba4ae4c9da4e58af2c7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Sat, 11 Jan 2025 12:04:42 +0000 Subject: [PATCH 13/35] Remove final references to India Def --- finat/__init__.py | 2 +- finat/element_factory.py | 4 ++-- finat/fiat_elements.py | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/finat/__init__.py b/finat/__init__.py index c02fc90d6..4183d2d3f 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -8,7 +8,7 @@ TrimmedSerendipityCurl, TrimmedSerendipityDiv, # noqa: F401 TrimmedSerendipityEdge, TrimmedSerendipityFace, # noqa: F401 Nedelec, NedelecSecondKind, RaviartThomas, Regge, # noqa: F401 - IndiaDefElement) # noqa: F401 + FuseElement) # noqa: F401 from .argyris import Argyris # noqa: F401 from .aw import ArnoldWinther, ArnoldWintherNC # noqa: F401 diff --git a/finat/element_factory.py b/finat/element_factory.py index 5adcca5f9..a1947536e 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -306,8 +306,8 @@ def convert_restrictedelement(element, **kwargs): @convert.register(finat.ufl.FuseElement) -def convert_india_def(element, **kwargs): - return finat.fiat_elements.IndiaDefElement(element.triple), set() +def convert_fuse_element(element, **kwargs): + return finat.fiat_elements.FuseElement(element.triple), set() hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval) diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 4ac4aa7ba..e2b2d5107 100644 --- a/finat/fiat_elements.py +++ b/finat/fiat_elements.py @@ -482,6 +482,6 @@ def __init__(self, cell, degree, variant=None): super().__init__(FIAT.NedelecSecondKind(cell, degree, variant=variant)) -class IndiaDefElement(FiatElement): +class FuseElement(FiatElement): def __init__(self, triple): - super(IndiaDefElement, self).__init__(triple.to_fiat()) + super(FuseElement, self).__init__(triple.to_fiat()) From 2e3ed6c11a36e8f7ebcc9f35bff2d89691c336d6 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 13 Jan 2025 15:00:07 +0000 Subject: [PATCH 14/35] Add test, small change to spatial dimension management --- FIAT/polynomial_set.py | 8 +++--- test/FIAT/unit/test_polynomial.py | 43 +++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 4 deletions(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 4f61968bd..f9b9e8f7e 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -202,18 +202,18 @@ def construct_new_coeffs(ref_el, A, B): if A.get_expansion_set().continuity != B.get_expansion_set().continuity: raise ValueError("Continuity of expansion sets does not match.") if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None: - higher = A if A.degree > B.degree else B - lower = B if A.degree > B.degree else A + higher = A if A.get_embedded_degree() > B.get_embedded_degree() else B + lower = B if A.get_embedded_degree() > B.get_embedded_degree() else A try: sd = lower.get_shape()[0] except IndexError: - sd = 1 + sd = -1 embedded_coeffs = [] diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1] for coeff in lower.coeffs: - if sd > 1: + if sd != -1: new_coeff = [] for row in coeff: new_coeff.append(numpy.append(row, [0 for i in range(diff)])) diff --git a/test/FIAT/unit/test_polynomial.py b/test/FIAT/unit/test_polynomial.py index 2982cf11e..c211fb815 100644 --- a/test/FIAT/unit/test_polynomial.py +++ b/test/FIAT/unit/test_polynomial.py @@ -21,6 +21,7 @@ from FIAT import expansions, polynomial_set, reference_element from FIAT.quadrature_schemes import create_quadrature +from itertools import chain @pytest.fixture(params=(1, 2, 3)) @@ -107,3 +108,45 @@ def test_bubble_duality(cell, degree): results = scale * numpy.dot(numpy.multiply(phi_dual, qwts), phi.T) assert numpy.allclose(results, numpy.diag(numpy.diag(results))) assert numpy.allclose(numpy.diag(results), 1.0) + + +@pytest.mark.parametrize("degree", [10]) +def test_union_of_polysets(cell, degree): + """ demonstrates that polysets don't need to have the same degree for union + using RT space as an example""" + + sd = cell.get_spatial_dimension() + k = degree + vecPk = polynomial_set.ONPolynomialSet(cell, degree, (sd,)) + + vec_Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, (sd,), scale="orthonormal") + + dimPkp1 = expansions.polynomial_dimension(cell, k + 1) + dimPk = expansions.polynomial_dimension(cell, k) + dimPkm1 = expansions.polynomial_dimension(cell, k - 1) + + vec_Pk_indices = list(chain(*(range(i * dimPkp1, i * dimPkp1 + dimPk) + for i in range(sd)))) + vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) + + Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, scale="orthonormal") + PkH = Pkp1.take(list(range(dimPkm1, dimPk))) + + Q = create_quadrature(cell, 2 * (k + 1)) + Qpts, Qwts = Q.get_points(), Q.get_weights() + + PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] + x = Qpts.T + PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] + PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) + PkHx = polynomial_set.PolynomialSet(cell, k, k + 1, vec_Pkp1.get_expansion_set(), PkHx_coeffs) + + same_deg = polynomial_set.polynomial_set_union_normalized(vec_Pk_from_Pkp1, PkHx) + different_deg = polynomial_set.polynomial_set_union_normalized(vecPk, PkHx) + + Q = create_quadrature(cell, 2*(degree)) + Qpts, _ = Q.get_points(), Q.get_weights() + same_vals = same_deg.tabulate(Qpts)[(0,) * sd] + diff_vals = different_deg.tabulate(Qpts)[(0,) * sd] + assert numpy.allclose(same_vals - diff_vals, 0) From 0b593e556b681ef4b02e252a39579683792c542d Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 15 Jan 2025 18:08:56 +0000 Subject: [PATCH 15/35] remove fuse dependency --- finat/element_factory.py | 7 ++++--- pyproject.toml | 1 - 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/finat/element_factory.py b/finat/element_factory.py index a1947536e..154336764 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -24,7 +24,6 @@ import finat import finat.ufl -from fuse.cells import CellComplexToUFL as FuseCell import ufl from FIAT import ufc_cell @@ -113,9 +112,11 @@ def as_fiat_cell(cell): :arg cell: the :class:`ufl.Cell` to convert.""" if not isinstance(cell, ufl.AbstractCell): raise ValueError("Expecting a UFL Cell") - if isinstance(cell, FuseCell): + + try: return cell.to_fiat() - return ufc_cell(cell) + except AttributeError: + return ufc_cell(cell) @singledispatch diff --git a/pyproject.toml b/pyproject.toml index 2d43b6dc5..c224c1fbe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,6 @@ dependencies = [ "symengine", "sympy", "fenics-ufl @ git+https://github.com/firedrakeproject/ufl.git", - "fuse-element @ git+https://github.com/indiamai/fuse.git", ] requires-python = ">=3.10" authors = [ From 263dad1326f2d56af4c37ebdfc6516a1a8705cc7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 16 Jan 2025 18:07:48 +0000 Subject: [PATCH 16/35] Modify code to use np.pad --- FIAT/polynomial_set.py | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index f9b9e8f7e..9686cf245 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -198,32 +198,24 @@ def construct_new_coeffs(ref_el, A, B): This does not handle the case that A and B have continuity but not the same degree. """ - sd = ref_el.get_spatial_dimension() if A.get_expansion_set().continuity != B.get_expansion_set().continuity: raise ValueError("Continuity of expansion sets does not match.") + if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None: higher = A if A.get_embedded_degree() > B.get_embedded_degree() else B lower = B if A.get_embedded_degree() > B.get_embedded_degree() else A - try: - sd = lower.get_shape()[0] - except IndexError: - sd = -1 - - embedded_coeffs = [] diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1] - for coeff in lower.coeffs: - if sd != -1: - new_coeff = [] - for row in coeff: - new_coeff.append(numpy.append(row, [0 for i in range(diff)])) - embedded_coeffs.append(new_coeff) - else: - embedded_coeffs.append(numpy.append(coeff, [0 for i in range(diff)])) - embedded_coeffs = numpy.array(embedded_coeffs) + + # pad only the 0th axis with the difference in size + padding = [(0, 0) for i in range(len(lower.coeffs.shape) - 1)] + [(0, diff)] + embedded_coeffs = numpy.pad(lower.coeffs, padding) + new_coeffs = numpy.concatenate((embedded_coeffs, higher.coeffs), axis=0) - else: + elif A.get_embedded_degree() == B.get_embedded_degree(): new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + else: + raise NotImplementedError("Extending of coefficients is not implemented for PolynomialSets with continuity and different degrees") return new_coeffs From ad63ee00ff4d2436394a0272455ac5ca53a0f042 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 16 Jan 2025 18:41:20 +0000 Subject: [PATCH 17/35] Add extension of coefficients to allow union of polysets with different degrees --- FIAT/nedelec.py | 8 +++--- FIAT/polynomial_set.py | 40 +++++++++++++++++++++++++--- FIAT/raviart_thomas.py | 6 ++--- test/FIAT/unit/test_polynomial.py | 43 +++++++++++++++++++++++++++++++ 4 files changed, 87 insertions(+), 10 deletions(-) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c92fa06be..7be2efcb8 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -22,7 +22,7 @@ def NedelecSpace2D(ref_el, degree): raise Exception("NedelecSpace2D requires 2d reference element") k = degree - 1 - vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,)) + vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal") dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1) dimPk = expansions.polynomial_dimension(ref_el, k) @@ -32,7 +32,7 @@ def NedelecSpace2D(ref_el, degree): for i in range(sd)))) vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) - Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1) + Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal") PkH = Pkp1.take(list(range(dimPkm1, dimPk))) Q = create_quadrature(ref_el, 2 * (k + 1)) @@ -43,8 +43,8 @@ def NedelecSpace2D(ref_el, degree): CrossX = numpy.dot(numpy.array([[0.0, 1.0], [-1.0, 0.0]]), Qpts.T) PkHCrossX_at_Qpts = PkH_at_Qpts[:, None, :] * CrossX[None, :, :] - PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), Pkp1_at_Qpts.T) - + PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), + Pkp1_at_Qpts.T) PkHcrossX = polynomial_set.PolynomialSet(ref_el, k + 1, k + 1, diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 08e62b8d8..9686cf245 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -176,15 +176,49 @@ def polynomial_set_union_normalized(A, B): not contain any of the same members of the set, as we construct a span via SVD. """ - new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + assert A.get_reference_element() == B.get_reference_element() + new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B) + + deg = max(A.get_degree(), B.get_degree()) + em_deg = max(A.get_embedded_degree(), B.get_embedded_degree()) coeffs = spanning_basis(new_coeffs) return PolynomialSet(A.get_reference_element(), - A.get_degree(), - A.get_embedded_degree(), + deg, + em_deg, A.get_expansion_set(), coeffs) +def construct_new_coeffs(ref_el, A, B): + """ + Constructs new coefficients for the set union of A and B + If A and B are discontinuous and do not have the same degree the smaller one + is extended to match the larger. + + This does not handle the case that A and B have continuity but not the same degree. + """ + + if A.get_expansion_set().continuity != B.get_expansion_set().continuity: + raise ValueError("Continuity of expansion sets does not match.") + + if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None: + higher = A if A.get_embedded_degree() > B.get_embedded_degree() else B + lower = B if A.get_embedded_degree() > B.get_embedded_degree() else A + + diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1] + + # pad only the 0th axis with the difference in size + padding = [(0, 0) for i in range(len(lower.coeffs.shape) - 1)] + [(0, diff)] + embedded_coeffs = numpy.pad(lower.coeffs, padding) + + new_coeffs = numpy.concatenate((embedded_coeffs, higher.coeffs), axis=0) + elif A.get_embedded_degree() == B.get_embedded_degree(): + new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + else: + raise NotImplementedError("Extending of coefficients is not implemented for PolynomialSets with continuity and different degrees") + return new_coeffs + + class ONSymTensorPolynomialSet(PolynomialSet): """Constructs an orthonormal basis for symmetric-tensor-valued polynomials on a reference element. diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index db0058d12..2b9241899 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -20,7 +20,7 @@ def RTSpace(ref_el, degree): sd = ref_el.get_spatial_dimension() k = degree - 1 - vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,)) + vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal") dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1) dimPk = expansions.polynomial_dimension(ref_el, k) @@ -30,7 +30,7 @@ def RTSpace(ref_el, degree): for i in range(sd)))) vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) - Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1) + Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal") PkH = Pkp1.take(list(range(dimPkm1, dimPk))) Q = create_quadrature(ref_el, 2 * (k + 1)) @@ -39,12 +39,12 @@ def RTSpace(ref_el, degree): # have to work on this through "tabulate" interface # first, tabulate PkH at quadrature points PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] x = Qpts.T PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) - PkHx = polynomial_set.PolynomialSet(ref_el, k, k + 1, diff --git a/test/FIAT/unit/test_polynomial.py b/test/FIAT/unit/test_polynomial.py index 2982cf11e..c211fb815 100644 --- a/test/FIAT/unit/test_polynomial.py +++ b/test/FIAT/unit/test_polynomial.py @@ -21,6 +21,7 @@ from FIAT import expansions, polynomial_set, reference_element from FIAT.quadrature_schemes import create_quadrature +from itertools import chain @pytest.fixture(params=(1, 2, 3)) @@ -107,3 +108,45 @@ def test_bubble_duality(cell, degree): results = scale * numpy.dot(numpy.multiply(phi_dual, qwts), phi.T) assert numpy.allclose(results, numpy.diag(numpy.diag(results))) assert numpy.allclose(numpy.diag(results), 1.0) + + +@pytest.mark.parametrize("degree", [10]) +def test_union_of_polysets(cell, degree): + """ demonstrates that polysets don't need to have the same degree for union + using RT space as an example""" + + sd = cell.get_spatial_dimension() + k = degree + vecPk = polynomial_set.ONPolynomialSet(cell, degree, (sd,)) + + vec_Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, (sd,), scale="orthonormal") + + dimPkp1 = expansions.polynomial_dimension(cell, k + 1) + dimPk = expansions.polynomial_dimension(cell, k) + dimPkm1 = expansions.polynomial_dimension(cell, k - 1) + + vec_Pk_indices = list(chain(*(range(i * dimPkp1, i * dimPkp1 + dimPk) + for i in range(sd)))) + vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) + + Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, scale="orthonormal") + PkH = Pkp1.take(list(range(dimPkm1, dimPk))) + + Q = create_quadrature(cell, 2 * (k + 1)) + Qpts, Qwts = Q.get_points(), Q.get_weights() + + PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] + x = Qpts.T + PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] + PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) + PkHx = polynomial_set.PolynomialSet(cell, k, k + 1, vec_Pkp1.get_expansion_set(), PkHx_coeffs) + + same_deg = polynomial_set.polynomial_set_union_normalized(vec_Pk_from_Pkp1, PkHx) + different_deg = polynomial_set.polynomial_set_union_normalized(vecPk, PkHx) + + Q = create_quadrature(cell, 2*(degree)) + Qpts, _ = Q.get_points(), Q.get_weights() + same_vals = same_deg.tabulate(Qpts)[(0,) * sd] + diff_vals = different_deg.tabulate(Qpts)[(0,) * sd] + assert numpy.allclose(same_vals - diff_vals, 0) From 3a5b073187c4540b75f54fa1ad5155ef12b7d135 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 17 Jan 2025 16:28:35 +0000 Subject: [PATCH 18/35] Small changes to allow different topologies --- FIAT/lagrange.py | 2 +- FIAT/reference_element.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index eead0cfa6..dab845446 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -40,7 +40,7 @@ def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=Fal entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])] if sort_entities: # sort the entities by support vertex ids - support = [top[dim][entity] for dim, entity in entities] + support = [sorted(top[dim][entity]) for dim, entity in entities] entities = [entity for verts, entity in sorted(zip(support, entities))] # make nodes by getting points diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 7d9ff2968..cbb7607e3 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -515,7 +515,8 @@ def make_points(self, dim, entity_id, order, variant=None): facet of dimension dim. Order indicates how many points to include in each direction.""" if dim == 0: - return (self.get_vertices()[entity_id], ) + return (self.get_vertices()[self.get_topology()[dim][entity_id][0]],) + # return (self.get_vertices()[entity_id], ) elif 0 < dim <= self.get_spatial_dimension(): entity_verts = \ self.get_vertices_of_subcomplex( From 8d52e227a3f5852d00d98bd1ff9146447f24b295 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 22 Jan 2025 14:19:58 +0000 Subject: [PATCH 19/35] remove orthonormal as seems to be no longer necessary --- FIAT/nedelec.py | 2 +- FIAT/raviart_thomas.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 7be2efcb8..db8b26b3b 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -22,7 +22,7 @@ def NedelecSpace2D(ref_el, degree): raise Exception("NedelecSpace2D requires 2d reference element") k = degree - 1 - vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal") + vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,)) dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1) dimPk = expansions.polynomial_dimension(ref_el, k) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 2b9241899..21b45d0b1 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -20,7 +20,7 @@ def RTSpace(ref_el, degree): sd = ref_el.get_spatial_dimension() k = degree - 1 - vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal") + vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,)) dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1) dimPk = expansions.polynomial_dimension(ref_el, k) From 8eaa60cf3ee656a4162cccbd2b99fb3697f5e045 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 23 Jan 2025 10:54:50 +0000 Subject: [PATCH 20/35] Remove further orthonormals --- FIAT/nedelec.py | 2 +- FIAT/raviart_thomas.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index db8b26b3b..6cf122772 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -32,7 +32,7 @@ def NedelecSpace2D(ref_el, degree): for i in range(sd)))) vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) - Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal") + Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1) PkH = Pkp1.take(list(range(dimPkm1, dimPk))) Q = create_quadrature(ref_el, 2 * (k + 1)) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 21b45d0b1..2c7d72c83 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -30,7 +30,7 @@ def RTSpace(ref_el, degree): for i in range(sd)))) vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) - Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal") + Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1) PkH = Pkp1.take(list(range(dimPkm1, dimPk))) Q = create_quadrature(ref_el, 2 * (k + 1)) From f31fd72e731f395f1dbb216a6d40482e445253c6 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 10 Feb 2025 17:42:02 +0000 Subject: [PATCH 21/35] Update UFL rep of fuse element to cover both regular fuse and tensor prod --- finat/ufl/fuseelement.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py index b83b26ba0..ca79d8494 100644 --- a/finat/ufl/fuseelement.py +++ b/finat/ufl/fuseelement.py @@ -19,13 +19,12 @@ def __init__(self, triple, cell=None): if not cell: cell = self.triple.cell.to_ufl() - # this isn't really correct - degree = self.triple.spaces[0].degree() + degree = self.triple.degree() + self.sobolev_space = self.triple.spaces[1] super(FuseElement, self).__init__("IT", cell, degree, None, triple.get_value_shape()) def __repr__(self): - return "FiniteElement(%s, %s, (%s, %s, %s), %s)" % ( - repr(self.triple.DOFGenerator), repr(self.triple.cell), repr(self.triple.spaces[0]), repr(self.triple.spaces[1]), repr(self.triple.spaces[2]), "X") + return repr(self.triple) def __str__(self): return "" % (self.triple.spaces[0], self.triple.cell) From a91c7a2ca152d12133778ba381b3306601abac96 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 17 Feb 2025 16:36:55 +0000 Subject: [PATCH 22/35] tensor prod infrastructure --- finat/element_factory.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/finat/element_factory.py b/finat/element_factory.py index bd1e6cc57..7c23d6555 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -27,6 +27,7 @@ import ufl from FIAT import ufc_cell +from FIAT.reference_element import TensorProductCell __all__ = ("as_fiat_cell", "create_base_element", "create_element", "supported_elements") @@ -112,7 +113,8 @@ def as_fiat_cell(cell): :arg cell: the :class:`ufl.Cell` to convert.""" if not isinstance(cell, ufl.AbstractCell): raise ValueError("Expecting a UFL Cell") - + if isinstance(cell, ufl.TensorProductCell) and any([hasattr(c, "to_fiat") for c in cell._cells]): + return TensorProductCell(*[c.to_fiat() for c in cell._cells]) try: return cell.to_fiat() except AttributeError: From 2bc3161c5500a06981475d7bfa4b3e9e24d63b60 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 20 Feb 2025 17:47:15 +0000 Subject: [PATCH 23/35] convert things to generic as_Cell and add flattening infra for fuse --- FIAT/reference_element.py | 15 ++++++++++----- finat/cube.py | 6 +++--- finat/element_factory.py | 6 ++++++ 3 files changed, 19 insertions(+), 8 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index c71452594..d354ebee7 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1735,17 +1735,21 @@ def tuple_sum(tree): def is_hypercube(cell): + res = False if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)): - return True + res = True elif isinstance(cell, TensorProductCell): - return reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) + res = reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) else: - return False - + res = False + res2 = len(cell.vertices()) == 2 ** cell.get_dimension() + assert res == res2 + return res def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" - flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()} + from finat import as_fiat_cell + flattened_cube = {2: as_fiat_cell("quadrilateral"), 3: as_fiat_cell("hexahedron")} if numpy.sum(ref_el.get_dimension()) <= 1: # Just return point/interval cell arguments return ref_el @@ -1798,6 +1802,7 @@ def compute_unflattening_map(topology_dict): def max_complex(complexes): + breakpoint() max_cell = max(complexes) if all(max_cell >= b for b in complexes): return max_cell diff --git a/finat/cube.py b/finat/cube.py index 39aeccd04..d529e4237 100644 --- a/finat/cube.py +++ b/finat/cube.py @@ -5,7 +5,7 @@ flatten_permutations) from FIAT.tensor_product import FlattenedDimensions as FIAT_FlattenedDimensions from gem.utils import cached_property - +from ufl import as_cell from finat.finiteelementbase import FiniteElementBase @@ -23,9 +23,9 @@ def __init__(self, element): def cell(self): dim = self.product.cell.get_spatial_dimension() if dim == 2: - return UFCQuadrilateral() + return as_fiat_cell("quadrilateral") elif dim == 3: - return UFCHexahedron() + return as_fiat_cell("hexahedron") else: raise NotImplementedError("Cannot guess cell for spatial dimension %s" % dim) diff --git a/finat/element_factory.py b/finat/element_factory.py index 7c23d6555..df030219a 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -113,6 +113,8 @@ def as_fiat_cell(cell): :arg cell: the :class:`ufl.Cell` to convert.""" if not isinstance(cell, ufl.AbstractCell): raise ValueError("Expecting a UFL Cell") + if isinstance(cell, str): + cell = ufl.as_cell(cell) if isinstance(cell, ufl.TensorProductCell) and any([hasattr(c, "to_fiat") for c in cell._cells]): return TensorProductCell(*[c.to_fiat() for c in cell._cells]) try: @@ -313,6 +315,10 @@ def convert_restrictedelement(element, **kwargs): @convert.register(finat.ufl.FuseElement) def convert_fuse_element(element, **kwargs): + if element.triple.flat: + new_elem = element.triple.unflatten() + finat_elem, deps = _create_element(new_elem.to_ufl(), **kwargs) + return finat.FlattenedDimensions(finat_elem), deps return finat.fiat_elements.FuseElement(element.triple), set() From 4953241f971de86af98d25eb51154b233c80f2c7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 21 Feb 2025 16:05:58 +0000 Subject: [PATCH 24/35] Add generic hypercube class --- FIAT/reference_element.py | 88 +++++++++++++++++++++++++++++++++++++++ finat/cube.py | 5 +-- finat/element_factory.py | 8 ++-- 3 files changed, 94 insertions(+), 7 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index d354ebee7..ea6a2a8cf 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1344,6 +1344,93 @@ def extrinsic_orientation_permutation_map(self): return a +class Hypercube(Cell): + """ + For reference elements based on TensorProductCells""" + def __init__(self, product): + pt = product.get_topology() + + verts = product.get_vertices() + topology = flatten_entities(pt) + + # TODO this should be generalised ? + cube_types = {2: QUADRILATERAL, 3: HEXAHEDRON} + super().__init__(cube_types[sum(product.get_dimension())], verts, topology) + + self.product = product + self.unflattening_map = compute_unflattening_map(pt) + + def get_dimension(self): + """Returns the subelement dimension of the cell. Same as the + spatial dimension.""" + return self.get_spatial_dimension() + + def get_entity_transform(self, dim, entity_i): + """Returns a mapping of point coordinates from the + `entity_i`-th subentity of dimension `dim` to the cell. + + :arg dim: entity dimension (integer) + :arg entity_i: entity number (integer) + """ + d, e = self.unflattening_map[(dim, entity_i)] + return self.product.get_entity_transform(d, e) + + def volume(self): + """Computes the volume in the appropriate dimensional measure.""" + return self.product.volume() + + def compute_reference_normal(self, facet_dim, facet_i): + """Returns the unit normal in infinity norm to facet_i.""" + assert facet_dim == 1 + d, i = self.unflattening_map[(facet_dim, facet_i)] + return self.product.compute_reference_normal(d, i) + + def contains_point(self, point, epsilon=0): + """Checks if reference cell contains given point + (with numerical tolerance as given by the L1 distance (aka 'manhatten', + 'taxicab' or rectilinear distance) to the cell. + + Parameters + ---------- + point : numpy.ndarray, list or symbolic expression + The coordinates of the point. + epsilon : float + The tolerance for the check. + + Returns + ------- + bool : True if the point is inside the cell, False otherwise. + + """ + return self.product.contains_point(point, epsilon=epsilon) + + def distance_to_point_l1(self, point, rescale=False): + """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear + distance) to a point with 0.0 if the point is inside the cell. + + For more information see the docstring for the UFCSimplex method.""" + return self.product.distance_to_point_l1(point, rescale=rescale) + + def symmetry_group_size(self, dim): + return [1, 2, 8][dim] + + def cell_orientation_reflection_map(self): + """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" + return self.product.cell_orientation_reflection_map() + + def __gt__(self, other): + return self.product > other + + def __lt__(self, other): + return self.product < other + + def __ge__(self, other): + return self.product >= other + + def __le__(self, other): + return self.product <= other + + class UFCQuadrilateral(Cell): r"""This is the reference quadrilateral with vertices (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). @@ -1746,6 +1833,7 @@ def is_hypercube(cell): assert res == res2 return res + def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" from finat import as_fiat_cell diff --git a/finat/cube.py b/finat/cube.py index d529e4237..45bf206a7 100644 --- a/finat/cube.py +++ b/finat/cube.py @@ -1,12 +1,11 @@ from __future__ import absolute_import, division, print_function -from FIAT.reference_element import (UFCHexahedron, UFCQuadrilateral, - compute_unflattening_map, flatten_entities, +from FIAT.reference_element import (compute_unflattening_map, flatten_entities, flatten_permutations) from FIAT.tensor_product import FlattenedDimensions as FIAT_FlattenedDimensions from gem.utils import cached_property -from ufl import as_cell from finat.finiteelementbase import FiniteElementBase +from finat.element_factory import as_fiat_cell class FlattenedDimensions(FiniteElementBase): diff --git a/finat/element_factory.py b/finat/element_factory.py index df030219a..92bff89ba 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -111,10 +111,10 @@ def as_fiat_cell(cell): """Convert a ufl cell to a FIAT cell. :arg cell: the :class:`ufl.Cell` to convert.""" - if not isinstance(cell, ufl.AbstractCell): - raise ValueError("Expecting a UFL Cell") if isinstance(cell, str): cell = ufl.as_cell(cell) + if not isinstance(cell, ufl.AbstractCell): + raise ValueError("Expecting a UFL Cell") if isinstance(cell, ufl.TensorProductCell) and any([hasattr(c, "to_fiat") for c in cell._cells]): return TensorProductCell(*[c.to_fiat() for c in cell._cells]) try: @@ -322,8 +322,8 @@ def convert_fuse_element(element, **kwargs): return finat.fiat_elements.FuseElement(element.triple), set() -hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval) -quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval) +hexahedron_tpc = ufl.TensorProductCell(ufl.as_cell("interval"), ufl.as_cell("interval"), ufl.as_cell("interval")) +quadrilateral_tpc = ufl.TensorProductCell(ufl.as_cell("interval"), ufl.as_cell("interval")) _cache = weakref.WeakKeyDictionary() From e8aaa56364333aa9269bf04a587f963d2a06c452 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 26 Feb 2025 11:30:07 +0000 Subject: [PATCH 25/35] Change type to allow comparision of subclasses of tensor product cell with tensorproduct cell --- FIAT/reference_element.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index ea6a2a8cf..4debb4a36 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1257,7 +1257,7 @@ def cell_orientation_reflection_map(self): def compare(self, op, other): if hasattr(other, "product"): other = other.product - if isinstance(other, type(self)): + if isinstance(other, TensorProductCell): return all(op(a, b) for a, b in zip(self.cells, other.cells)) else: return op(self, other) @@ -1890,7 +1890,6 @@ def compute_unflattening_map(topology_dict): def max_complex(complexes): - breakpoint() max_cell = max(complexes) if all(max_cell >= b for b in complexes): return max_cell From a726c98b16e88e636f8402d8b6d40614c01048a1 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 13 Mar 2025 17:45:05 +0000 Subject: [PATCH 26/35] Add function that moves any cell to a simplex. Modify DPC to use --- FIAT/discontinuous_pc.py | 20 +++++++------------- FIAT/reference_element.py | 16 ++++++++++++---- 2 files changed, 19 insertions(+), 17 deletions(-) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index dbb566906..ba9894c91 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -15,21 +15,15 @@ UFCTriangle, UFCTetrahedron, make_affine_mapping, - flatten_reference_cube) + flatten_reference_cube, + cell_to_simplex) from FIAT.P0 import P0Dual import numpy as np -hypercube_simplex_map = {Point(): Point(), - DefaultLine(): DefaultLine(), - UFCInterval(): UFCInterval(), - UFCQuadrilateral(): UFCTriangle(), - UFCHexahedron(): UFCTetrahedron()} - - class DPC0(finite_element.CiarletElement): def __init__(self, ref_el): flat_el = flatten_reference_cube(ref_el) - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], 0) + poly_set = polynomial_set.ONPolynomialSet(cell_to_simplex(flat_el), 0) dual = P0Dual(ref_el) # Implement entity_permutations when we handle that for HigherOrderDPC. # Currently, orientation_tuples in P0Dual(ref_el).entity_permutations @@ -58,7 +52,7 @@ def __init__(self, ref_el, flat_el, degree): # Change coordinates here. # Vertices of the simplex corresponding to the reference element. - v_simplex = hypercube_simplex_map[flat_el].get_vertices() + v_simplex = cell_to_simplex(flat_el).get_vertices() # Vertices of the reference element. v_hypercube = flat_el.get_vertices() # For the mapping, first two vertices are unchanged in all dimensions. @@ -74,12 +68,12 @@ def __init__(self, ref_el, flat_el, degree): # make nodes by getting points # need to do this dimension-by-dimension, facet-by-facet - top = hypercube_simplex_map[flat_el].get_topology() + top = cell_to_simplex(flat_el).get_topology() cur = 0 for dim in sorted(top): for entity in sorted(top[dim]): - pts_cur = hypercube_simplex_map[flat_el].make_points(dim, entity, degree) + pts_cur = cell_to_simplex(flat_el).make_points(dim, entity, degree) pts_cur = [tuple(np.matmul(A, np.array(x)) + b) for x in pts_cur] nodes_cur = [functional.PointEvaluation(flat_el, x) for x in pts_cur] @@ -102,7 +96,7 @@ class HigherOrderDPC(finite_element.CiarletElement): def __init__(self, ref_el, degree): flat_el = flatten_reference_cube(ref_el) - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], degree) + poly_set = polynomial_set.ONPolynomialSet(cell_to_simplex(flat_el), degree) dual = DPCDualSet(ref_el, flat_el, degree) formdegree = flat_el.get_spatial_dimension() # n-form super().__init__(poly_set=poly_set, diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 4debb4a36..606275656 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1829,14 +1829,15 @@ def is_hypercube(cell): res = reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) else: res = False - res2 = len(cell.vertices()) == 2 ** cell.get_dimension() - assert res == res2 - return res + # breakpoint() + res2 = len(cell.vertices) == 2 ** cell.get_dimension() + # assert res == res2 + return res2 def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" - from finat import as_fiat_cell + from finat.element_factory import as_fiat_cell flattened_cube = {2: as_fiat_cell("quadrilateral"), 3: as_fiat_cell("hexahedron")} if numpy.sum(ref_el.get_dimension()) <= 1: # Just return point/interval cell arguments @@ -1895,3 +1896,10 @@ def max_complex(complexes): return max_cell else: raise ValueError("Cannot find the maximal complex") + + +def cell_to_simplex(cell): + if cell.is_simplex(): + return cell + else: + return ufc_simplex(cell.get_dimension()) From 7d5ef1c848446dba25b3ab4a50726e3a892a5d1d Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Fri, 21 Mar 2025 11:07:53 +0000 Subject: [PATCH 27/35] Integrate hypercube changes to fuse (#137) Co-authored-by: Pablo Brubeck Co-authored-by: Connor Ward --- FIAT/discontinuous_pc.py | 10 +- FIAT/dual_set.py | 3 + FIAT/expansions.py | 34 ++- FIAT/finite_element.py | 3 + FIAT/macro.py | 59 ++-- FIAT/polynomial_set.py | 3 + FIAT/reference_element.py | 334 ++++++++--------------- finat/point_set.py | 35 ++- finat/quadrature.py | 36 ++- gem/utils.py | 55 ++++ test/FIAT/unit/test_macro.py | 21 ++ test/FIAT/unit/test_reference_element.py | 59 ++++ test/finat/test_create_fiat_element.py | 5 + test/finat/test_quadrature.py | 19 ++ 14 files changed, 401 insertions(+), 275 deletions(-) create mode 100644 test/finat/test_quadrature.py diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index ba9894c91..8973a58f6 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -7,19 +7,13 @@ # Modified by David A. Ham (david.ham@imperial.ac.uk), 2018 from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.reference_element import (Point, - DefaultLine, - UFCInterval, - UFCQuadrilateral, - UFCHexahedron, - UFCTriangle, - UFCTetrahedron, - make_affine_mapping, +from FIAT.reference_element import (make_affine_mapping, flatten_reference_cube, cell_to_simplex) from FIAT.P0 import P0Dual import numpy as np + class DPC0(finite_element.CiarletElement): def __init__(self, ref_el): flat_el = flatten_reference_cube(ref_el) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 74630363d..c09db974f 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -35,6 +35,9 @@ def __init__(self, nodes, ref_el, entity_ids, entity_permutations=None): def __iter__(self): return iter(self.nodes) + def __len__(self): + return len(self.nodes) + def get_nodes(self): return self.nodes diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 28eb561d4..cf0a67a3b 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -436,26 +436,38 @@ def tabulate_jumps(self, n, points, order=0): """ from FIAT.polynomial_set import mis + sd = self.ref_el.get_spatial_dimension() num_members = self.get_num_members(n) cell_node_map = self.get_cell_node_map(n) + cell_point_map = compute_cell_point_map(self.ref_el, points, unique=False) + + num_jumps = 0 + facet_point_map = {} + for facet in self.ref_el.get_interior_facets(sd-1): + try: + cells = self.ref_el.connectivity[(sd-1, sd)][facet] + ipts = list(set.intersection(*(set(cell_point_map[c]) for c in cells))) + if ipts != (): + facet_point_map[facet] = ipts + num_jumps += len(ipts) + except KeyError: + pass - top = self.ref_el.get_topology() - sd = self.ref_el.get_spatial_dimension() - interior_facets = self.ref_el.get_interior_facets(sd-1) derivs = {cell: self._tabulate_on_cell(n, points, order=order, cell=cell) - for cell in top[sd]} + for cell in cell_point_map} + jumps = {} for r in range(order+1): cur = 0 alphas = mis(sd, r) - jumps[r] = numpy.zeros((num_members, len(alphas)*len(interior_facets), len(points))) - for facet in interior_facets: - facet_verts = set(top[sd-1][facet]) - c0, c1 = tuple(k for k in top[sd] if facet_verts < set(top[sd][k])) + jumps[r] = numpy.zeros((num_members, len(alphas) * num_jumps)) + for facet, ipts in facet_point_map.items(): + c0, c1 = self.ref_el.connectivity[(sd-1, sd)][facet] for alpha in alphas: - jumps[r][cell_node_map[c1], cur] += derivs[c1][alpha] - jumps[r][cell_node_map[c0], cur] -= derivs[c0][alpha] - cur = cur + 1 + ijump = range(cur, cur + len(ipts)) + jumps[r][numpy.ix_(cell_node_map[c1], ijump)] += derivs[c1][alpha][:, ipts] + jumps[r][numpy.ix_(cell_node_map[c0], ijump)] -= derivs[c0][alpha][:, ipts] + cur += len(ipts) return jumps def get_dmats(self, degree, cell=0): diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 6f4ff24d5..96cf26286 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -134,6 +134,9 @@ def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref ref_complex = ref_complex or poly_set.get_reference_element() super().__init__(ref_el, dual, order, formdegree, mapping, ref_complex) + if len(poly_set) != len(dual): + raise ValueError(f"Dimension of function space is {len(poly_set)}, but got {len(dual)} nodes.") + # build generalized Vandermonde matrix old_coeffs = poly_set.get_coeffs() dualmat = dual.to_riesz(poly_set) diff --git a/FIAT/macro.py b/FIAT/macro.py index 6662bdfcb..07a809acb 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -412,21 +412,35 @@ class CkPolynomialSet(polynomial_set.PolynomialSet): :arg ref_el: The simplicial complex. :arg degree: The polynomial degree. - :kwarg order: The order of continuity across subcells. + :kwarg order: The order of continuity across facets or a dict of dicts + mapping dimension and entity_id to the order of continuity at each entity. :kwarg vorder: The order of super-smoothness at interior vertices. :kwarg shape: The value shape. :kwarg variant: The variant for the underlying ExpansionSet. :kwarg scale: The scale for the underlying ExpansionSet. """ - def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): + def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): from FIAT.quadrature_schemes import create_quadrature + + if not isinstance(order, (int, dict)): + raise TypeError(f"'order' must be either an int or dict, not {type(order).__name__}") + + sd = ref_el.get_spatial_dimension() + if isinstance(order, int): + order = {sd-1: dict.fromkeys(ref_el.get_interior_facets(sd-1), order)} + if vorder is not None: + order[0] = dict.fromkeys(ref_el.get_interior_facets(0), vorder) + elif 0 not in order: + order[0] = {} + + if not all(k in {0, sd-1} for k in order): + raise NotImplementedError("Only face or vertex constraints have been implemented.") + expansion_set = expansions.ExpansionSet(ref_el, **kwargs) k = 1 if expansion_set.continuity == "C0" else 0 - num_members = expansion_set.get_num_members(degree) - sd = ref_el.get_spatial_dimension() + # Impose C^forder continuity across interior facets facet_el = ref_el.construct_subelement(sd-1) - phi_deg = 0 if sd == 1 else degree - k phi = polynomial_set.ONPolynomialSet(facet_el, phi_deg) Q = create_quadrature(facet_el, 2 * phi_deg) @@ -435,30 +449,33 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): weights = numpy.multiply(phi_at_qpts, qwts) rows = [] - for facet in ref_el.get_interior_facets(sd-1): - jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order) - for r in range(k, order+1): + for facet in order[sd-1]: + forder = order[sd-1][facet] + jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=forder) + for r in range(k, forder+1): num_wt = 1 if sd == 1 else expansions.polynomial_dimension(facet_el, degree-r) - rows.append(numpy.tensordot(weights[:num_wt], jumps[r], axes=(-1, -1)).reshape(-1, num_members)) + rows.append(numpy.tensordot(weights[:num_wt], jumps[r], axes=(-1, -1))) + + # Impose C^vorder super-smoothness at interior vertices + # C^forder automatically gives C^{forder+dim-1} at the interior vertex + verts = numpy.asarray(ref_el.get_vertices()) + for vorder in set(order[0].values()): + vids = [i for i in order[0] if order[0][i] == vorder] + facets = chain.from_iterable(ref_el.connectivity[(0, sd-1)][v] for v in vids) + forder = min(order[sd-1][f] for f in facets) + sorder = forder + sd - 1 + if vorder > sorder: + jumps = expansion_set.tabulate_jumps(degree, verts[vids], order=vorder) + rows.extend(numpy.vstack(jumps[r].T) for r in range(sorder+1, vorder+1)) if len(rows) > 0: + for row in rows: + row *= 1 / max(numpy.max(abs(row)), 1) dual_mat = numpy.vstack(rows) coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True) else: coeffs = numpy.eye(expansion_set.get_num_members(degree)) - # Impose C^vorder super-smoothness at interior vertices - # C^order automatically gives C^{order+dim-1} at the interior vertex - sorder = order + sd - 1 - if vorder > sorder: - verts = ref_el.get_vertices() - points = [verts[i] for i in ref_el.get_interior_facets(0)] - jumps = expansion_set.tabulate_jumps(degree, points, order=vorder) - for r in range(sorder+1, vorder+1): - dual_mat = numpy.dot(numpy.vstack(jumps[r].T), coeffs.T) - nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) - coeffs = numpy.dot(nsp, coeffs) - if shape != (): m, n = coeffs.shape ncomp = numpy.prod(shape) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 9686cf245..534c18642 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -104,6 +104,9 @@ def take(self, items): return PolynomialSet(self.ref_el, self.degree, self.embedded_degree, self.expansion_set, new_coeffs) + def __len__(self): + return self.num_members + class ONPolynomialSet(PolynomialSet): """Constructs an orthonormal basis out of expansion set by having an diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 606275656..edc902f0d 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -26,6 +26,7 @@ from math import factorial import numpy +from gem.utils import safe_repr from recursivenodes.nodes import _decode_family, _recursive from FIAT.orientation_utils import ( @@ -43,6 +44,8 @@ HEXAHEDRON = 111 TENSORPRODUCT = 99 +hypercube_shapes = {2: QUADRILATERAL, 3: HEXAHEDRON} + def multiindex_equal(d, isum, imin=0): """A generator for d-tuple multi-indices whose sum is isum and minimum is imin. @@ -126,7 +129,7 @@ def linalg_subspace_intersection(A, B): return U[:, :rank_c] -class Cell(object): +class Cell: """Abstract class for a reference cell. Provides accessors for geometry (vertex coordinates) as well as topology (orderings of vertices that make up edges, faces, etc.""" @@ -184,6 +187,9 @@ def __init__(self, shape, vertices, topology): # Dictionary with derived cells self._split_cache = {} + def __repr__(self): + return f"{type(self).__name__}({self.shape!r}, {safe_repr(self.vertices)}, {self.topology!r})" + def _key(self): """Hashable object key data (excluding type).""" # Default: only type matters @@ -1131,6 +1137,9 @@ def __init__(self, *cells): super().__init__(TENSORPRODUCT, vertices, topology) self.cells = tuple(cells) + def __repr__(self): + return f"{type(self).__name__}({self.cells!r})" + def _key(self): return self.cells @@ -1212,7 +1221,7 @@ def compute_reference_normal(self, facet_dim, facet_i): def contains_point(self, point, epsilon=0.0): """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', + (with numerical tolerance as given by the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to the cell. Parameters @@ -1255,6 +1264,8 @@ def cell_orientation_reflection_map(self): return make_cell_orientation_reflection_map_tensorproduct(self.cells) def compare(self, op, other): + """Parent-based comparison between simplicial complexes. + This is done dimension by dimension.""" if hasattr(other, "product"): other = other.product if isinstance(other, TensorProductCell): @@ -1345,17 +1356,17 @@ def extrinsic_orientation_permutation_map(self): class Hypercube(Cell): - """ - For reference elements based on TensorProductCells""" - def __init__(self, product): - pt = product.get_topology() + """Abstract class for a reference hypercube""" + def __init__(self, dimension, product): + self.dimension = dimension + self.shape = hypercube_shapes[dimension] + + pt = product.get_topology() verts = product.get_vertices() topology = flatten_entities(pt) - # TODO this should be generalised ? - cube_types = {2: QUADRILATERAL, 3: HEXAHEDRON} - super().__init__(cube_types[sum(product.get_dimension())], verts, topology) + super().__init__(self.shape, verts, topology) self.product = product self.unflattening_map = compute_unflattening_map(pt) @@ -1365,6 +1376,21 @@ def get_dimension(self): spatial dimension.""" return self.get_spatial_dimension() + def construct_subelement(self, dimension): + """Constructs the reference element of a cell subentity + specified by subelement dimension. + + :arg dimension: subentity dimension (integer) + """ + sd = self.get_spatial_dimension() + if dimension > sd: + raise ValueError(f"Invalid dimension: {(dimension,)}") + elif dimension == sd: + return self + else: + sub_element = self.product.construct_subelement((dimension,) + (0,)*(len(self.product.cells) - 1)) + return flatten_reference_cube(sub_element) + def get_entity_transform(self, dim, entity_i): """Returns a mapping of point coordinates from the `entity_i`-th subentity of dimension `dim` to the cell. @@ -1381,13 +1407,14 @@ def volume(self): def compute_reference_normal(self, facet_dim, facet_i): """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 1 + sd = self.get_spatial_dimension() + assert facet_dim == sd - 1 d, i = self.unflattening_map[(facet_dim, facet_i)] return self.product.compute_reference_normal(d, i) def contains_point(self, point, epsilon=0): """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', + (with numerical tolerance as given by the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to the cell. Parameters @@ -1405,14 +1432,15 @@ def contains_point(self, point, epsilon=0): return self.product.contains_point(point, epsilon=epsilon) def distance_to_point_l1(self, point, rescale=False): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear + """Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. For more information see the docstring for the UFCSimplex method.""" return self.product.distance_to_point_l1(point, rescale=rescale) def symmetry_group_size(self, dim): - return [1, 2, 8][dim] + """Size of hypercube symmetry group is d! * 2**d""" + return factorial(dim) * (2**dim) def cell_orientation_reflection_map(self): """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" @@ -1431,7 +1459,33 @@ def __le__(self, other): return self.product <= other -class UFCQuadrilateral(Cell): +class UFCHypercube(Hypercube): + """Reference UFC Hypercube + + UFCHypercube: [0, 1]^d with vertices in + lexicographical order.""" + + def __init__(self, dim): + cells = [UFCInterval()] * dim + product = TensorProductCell(*cells) + super().__init__(dim, product) + + def construct_subelement(self, dimension): + """Constructs the reference element of a cell subentity + specified by subelement dimension. + + :arg dimension: subentity dimension (integer) + """ + sd = self.get_spatial_dimension() + if dimension > sd: + raise ValueError(f"Invalid dimension: {dimension}") + elif dimension == sd: + return self + else: + return ufc_hypercube(dimension) + + +class UFCQuadrilateral(UFCHypercube): r"""This is the reference quadrilateral with vertices (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). @@ -1476,207 +1530,18 @@ class UFCQuadrilateral(Cell): o = 2 * 1 + 0 = 2 """ - def __init__(self): - product = TensorProductCell(UFCInterval(), UFCInterval()) - pt = product.get_topology() - - verts = product.get_vertices() - topology = flatten_entities(pt) - - super().__init__(QUADRILATERAL, verts, topology) - - self.product = product - self.unflattening_map = compute_unflattening_map(pt) - - def get_dimension(self): - """Returns the subelement dimension of the cell. Same as the - spatial dimension.""" - return self.get_spatial_dimension() - - def construct_subelement(self, dimension): - """Constructs the reference element of a cell subentity - specified by subelement dimension. - - :arg dimension: subentity dimension (integer) - """ - if dimension == 2: - return self - elif dimension == 1: - return UFCInterval() - elif dimension == 0: - return Point() - else: - raise ValueError("Invalid dimension: %d" % (dimension,)) - - def get_entity_transform(self, dim, entity_i): - """Returns a mapping of point coordinates from the - `entity_i`-th subentity of dimension `dim` to the cell. - - :arg dim: entity dimension (integer) - :arg entity_i: entity number (integer) - """ - d, e = self.unflattening_map[(dim, entity_i)] - return self.product.get_entity_transform(d, e) - - def volume(self): - """Computes the volume in the appropriate dimensional measure.""" - return self.product.volume() - - def compute_reference_normal(self, facet_dim, facet_i): - """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 1 - d, i = self.unflattening_map[(facet_dim, facet_i)] - return self.product.compute_reference_normal(d, i) - - def contains_point(self, point, epsilon=0): - """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', - 'taxicab' or rectilinear distance) to the cell. - - Parameters - ---------- - point : numpy.ndarray, list or symbolic expression - The coordinates of the point. - epsilon : float - The tolerance for the check. - - Returns - ------- - bool : True if the point is inside the cell, False otherwise. - - """ - return self.product.contains_point(point, epsilon=epsilon) - - def distance_to_point_l1(self, point, rescale=False): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear - distance) to a point with 0.0 if the point is inside the cell. - - For more information see the docstring for the UFCSimplex method.""" - return self.product.distance_to_point_l1(point, rescale=rescale) - - def symmetry_group_size(self, dim): - return [1, 2, 8][dim] - - def cell_orientation_reflection_map(self): - """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" - return self.product.cell_orientation_reflection_map() - - def __gt__(self, other): - return self.product > other - - def __lt__(self, other): - return self.product < other - - def __ge__(self, other): - return self.product >= other - def __le__(self, other): - return self.product <= other + def __init__(self): + super().__init__(2) -class UFCHexahedron(Cell): +class UFCHexahedron(UFCHypercube): """This is the reference hexahedron with vertices (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" def __init__(self): - product = TensorProductCell(UFCInterval(), UFCInterval(), UFCInterval()) - pt = product.get_topology() - - verts = product.get_vertices() - topology = flatten_entities(pt) - - super().__init__(HEXAHEDRON, verts, topology) - - self.product = product - self.unflattening_map = compute_unflattening_map(pt) - - def get_dimension(self): - """Returns the subelement dimension of the cell. Same as the - spatial dimension.""" - return self.get_spatial_dimension() - - def construct_subelement(self, dimension): - """Constructs the reference element of a cell subentity - specified by subelement dimension. - - :arg dimension: subentity dimension (integer) - """ - if dimension == 3: - return self - elif dimension == 2: - return UFCQuadrilateral() - elif dimension == 1: - return UFCInterval() - elif dimension == 0: - return Point() - else: - raise ValueError("Invalid dimension: %d" % (dimension,)) - - def get_entity_transform(self, dim, entity_i): - """Returns a mapping of point coordinates from the - `entity_i`-th subentity of dimension `dim` to the cell. - - :arg dim: entity dimension (integer) - :arg entity_i: entity number (integer) - """ - d, e = self.unflattening_map[(dim, entity_i)] - return self.product.get_entity_transform(d, e) - - def volume(self): - """Computes the volume in the appropriate dimensional measure.""" - return self.product.volume() - - def compute_reference_normal(self, facet_dim, facet_i): - """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 2 - d, i = self.unflattening_map[(facet_dim, facet_i)] - return self.product.compute_reference_normal(d, i) - - def contains_point(self, point, epsilon=0): - """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', - 'taxicab' or rectilinear distance) to the cell. - - Parameters - ---------- - point : numpy.ndarray, list or symbolic expression - The coordinates of the point. - epsilon : float - The tolerance for the check. - - Returns - ------- - bool : True if the point is inside the cell, False otherwise. - - """ - return self.product.contains_point(point, epsilon=epsilon) - - def distance_to_point_l1(self, point, rescale=False): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear - distance) to a point with 0.0 if the point is inside the cell. - - For more information see the docstring for the UFCSimplex method.""" - return self.product.distance_to_point_l1(point, rescale=rescale) - - def symmetry_group_size(self, dim): - return [1, 2, 8, 48][dim] - - def cell_orientation_reflection_map(self): - """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" - return self.product.cell_orientation_reflection_map() - - def __gt__(self, other): - return self.product > other - - def __lt__(self, other): - return self.product < other - - def __ge__(self, other): - return self.product >= other - - def __le__(self, other): - return self.product <= other + super().__init__(3) def make_affine_mapping(xs, ys): @@ -1715,6 +1580,21 @@ def make_affine_mapping(xs, ys): return A, b +def ufc_hypercube(spatial_dim): + """Factory function that maps spatial dimension to an instance of + the UFC reference hypercube of that dimension.""" + if spatial_dim == 0: + return Point() + elif spatial_dim == 1: + return UFCInterval() + elif spatial_dim == 2: + return UFCQuadrilateral() + elif spatial_dim == 3: + return UFCHexahedron() + else: + raise RuntimeError(f"Can't create UFC hypercube of dimension {spatial_dim}.") + + def default_simplex(spatial_dim): """Factory function that maps spatial dimension to an instance of the default reference simplex of that dimension.""" @@ -1727,7 +1607,7 @@ def default_simplex(spatial_dim): elif spatial_dim == 3: return DefaultTetrahedron() else: - raise RuntimeError("Can't create default simplex of dimension %s." % str(spatial_dim)) + raise RuntimeError(f"Can't create default simplex of dimension {spatial_dim}.") def ufc_simplex(spatial_dim): @@ -1742,7 +1622,7 @@ def ufc_simplex(spatial_dim): elif spatial_dim == 3: return UFCTetrahedron() else: - raise RuntimeError("Can't create UFC simplex of dimension %s." % str(spatial_dim)) + raise RuntimeError(f"Can't create UFC simplex of dimension {spatial_dim}.") def symmetric_simplex(spatial_dim): @@ -1782,7 +1662,7 @@ def ufc_cell(cell): elif celltype == "tetrahedron": return ufc_simplex(3) else: - raise RuntimeError("Don't know how to create UFC cell of type %s" % str(celltype)) + raise RuntimeError(f"Don't know how to create UFC cell of type {str(celltype)}") def volume(verts): @@ -1821,32 +1701,38 @@ def tuple_sum(tree): return tree +def is_ufc(cell): + if isinstance(cell, (Point, UFCInterval, UFCHypercube, UFCSimplex)): + return True + elif isinstance(cell, TensorProductCell): + return all(is_ufc(c) for c in cell.cells) + else: + return False + + def is_hypercube(cell): - res = False - if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)): - res = True + if isinstance(cell, (DefaultLine, UFCInterval, Hypercube)): + return True elif isinstance(cell, TensorProductCell): - res = reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) + return all(is_hypercube(c) for c in cell.cells) else: - res = False - # breakpoint() - res2 = len(cell.vertices) == 2 ** cell.get_dimension() - # assert res == res2 - return res2 + return False def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" - from finat.element_factory import as_fiat_cell - flattened_cube = {2: as_fiat_cell("quadrilateral"), 3: as_fiat_cell("hexahedron")} - if numpy.sum(ref_el.get_dimension()) <= 1: + if ref_el.get_spatial_dimension() <= 1: # Just return point/interval cell arguments return ref_el else: # Handle cases where cell is a quad/cube constructed from a tensor product or # an already flattened element - if is_hypercube(ref_el): - return flattened_cube[numpy.sum(ref_el.get_dimension())] + if isinstance(ref_el, TensorProductCell): + if is_ufc(ref_el): + return ufc_hypercube(ref_el.get_spatial_dimension()) + return Hypercube(ref_el.get_spatial_dimension(), ref_el) + elif is_hypercube(ref_el): + return ref_el else: raise TypeError('Invalid cell type') diff --git a/finat/point_set.py b/finat/point_set.py index 1497308c7..2a437c9c1 100644 --- a/finat/point_set.py +++ b/finat/point_set.py @@ -1,13 +1,15 @@ -from abc import ABCMeta, abstractproperty +import abc +import hashlib +from functools import cached_property from itertools import chain, product import numpy import gem -from gem.utils import cached_property +from gem.utils import safe_repr -class AbstractPointSet(metaclass=ABCMeta): +class AbstractPointSet(abc.ABC): """A way of specifying a known set of points, perhaps with some (tensor) structure. @@ -15,8 +17,15 @@ class AbstractPointSet(metaclass=ABCMeta): where point_set_shape is () for scalar, (N,) for N element vector, (N, M) for N x M matrix etc. """ + def __hash__(self): + return int.from_bytes(hashlib.md5(repr(self).encode()).digest(), byteorder="big") - @abstractproperty + @abc.abstractmethod + def __repr__(self): + pass + + @property + @abc.abstractmethod def points(self): """A flattened numpy array of points or ``UnknownPointsArray`` object with shape (# of points, point dimension).""" @@ -27,12 +36,14 @@ def dimension(self): _, dim = self.points.shape return dim - @abstractproperty + @property + @abc.abstractmethod def indices(self): """GEM indices with matching shape and extent to the structure of the point set.""" - @abstractproperty + @property + @abc.abstractmethod def expression(self): """GEM expression describing the points, with free indices ``self.indices`` and shape (point dimension,).""" @@ -53,6 +64,9 @@ def __init__(self, point): assert len(point.shape) == 1 self.point = point + def __repr__(self): + return f"{type(self).__name__}({safe_repr(self.point)})" + @cached_property def points(self): # Make sure we conform to the expected (# of points, point dimension) @@ -106,6 +120,9 @@ def __init__(self, points_expr): assert len(points_expr.shape) == 2 self._points_expr = points_expr + def __repr__(self): + return f"{type(self).__name__}({self._points_expr!r})" + @cached_property def points(self): return UnknownPointsArray(self._points_expr.shape) @@ -133,6 +150,9 @@ def __init__(self, points): assert len(points.shape) == 2 self.points = points + def __repr__(self): + return f"{type(self).__name__}({self.points!r})" + @cached_property def points(self): pass # set at initialisation @@ -177,6 +197,9 @@ class TensorPointSet(AbstractPointSet): def __init__(self, factors): self.factors = tuple(factors) + def __repr__(self): + return f"{type(self).__name__}({self.factors!r})" + @cached_property def points(self): return numpy.array([list(chain(*pt_tuple)) diff --git a/finat/quadrature.py b/finat/quadrature.py index ec6def127..10e4048b9 100644 --- a/finat/quadrature.py +++ b/finat/quadrature.py @@ -1,12 +1,13 @@ -from abc import ABCMeta, abstractproperty -from functools import reduce +import hashlib +from abc import ABCMeta, abstractmethod +from functools import cached_property, reduce import gem import numpy from FIAT.quadrature import GaussLegendreQuadratureLineRule from FIAT.quadrature_schemes import create_quadrature as fiat_scheme from FIAT.reference_element import LINE, QUADRILATERAL, TENSORPRODUCT -from gem.utils import cached_property +from gem.utils import safe_repr from finat.point_set import GaussLegendrePointSet, PointSet, TensorPointSet @@ -60,11 +61,23 @@ class AbstractQuadratureRule(metaclass=ABCMeta): """Abstract class representing a quadrature rule as point set and a corresponding set of weights.""" - @abstractproperty + def __hash__(self): + return int.from_bytes(hashlib.md5(repr(self).encode()).digest(), byteorder="big") + + def __eq__(self, other): + return type(other) is type(self) and repr(other) == repr(self) + + @abstractmethod + def __repr__(self): + pass + + @property + @abstractmethod def point_set(self): """Point set object representing the quadrature points.""" - @abstractproperty + @property + @abstractmethod def weight_expression(self): """GEM expression describing the weights, with the same free indices as the point set.""" @@ -110,6 +123,16 @@ def __init__(self, point_set, weights, ref_el=None, io_ornt_map_tuple=(None, )): self.weights = numpy.asarray(weights) self._intrinsic_orientation_permutation_map_tuple = io_ornt_map_tuple + def __repr__(self): + return ( + f"{type(self).__name__}(" + f"{self.point_set!r}, " + f"{safe_repr(self.weights)}, " + f"{self.ref_el!r}, " + f"{self._intrinsic_orientation_permutation_map_tuple!r}" + ")" + ) + @cached_property def point_set(self): pass # set at initialisation @@ -131,6 +154,9 @@ def __init__(self, factors, ref_el=None): for m in factor._intrinsic_orientation_permutation_map_tuple ) + def __repr__(self): + return f"{type(self).__name__}({self.factors!r}, {self.ref_el!r})" + @cached_property def point_set(self): return TensorPointSet(q.point_set for q in self.factors) diff --git a/gem/utils.py b/gem/utils.py index 7e855f768..a6bf97dfa 100644 --- a/gem/utils.py +++ b/gem/utils.py @@ -1,5 +1,10 @@ import collections +import functools +import numbers from functools import cached_property # noqa: F401 +from typing import Any + +import numpy as np def groupby(iterable, key=None): @@ -88,3 +93,53 @@ def __exit__(self, exc_type, exc_value, traceback): assert self.state is variable._head value, variable._head = variable._head self.state = None + + +@functools.singledispatch +def safe_repr(obj: Any) -> str: + """Return a 'safe' repr for an object, accounting for floating point error. + + Parameters + ---------- + obj : + The object to produce a repr for. + + Returns + ------- + str : + A repr for the object. + + """ + raise TypeError(f"Cannot provide a safe repr for {type(obj).__name__}") + + +@safe_repr.register(str) +def _(text: str) -> str: + return text + + +@safe_repr.register(numbers.Integral) +def _(num: numbers.Integral) -> str: + return repr(num) + + +@safe_repr.register(numbers.Real) +def _(num: numbers.Real) -> str: + # set roundoff to close-to-but-not-exactly machine epsilon + precision = np.finfo(num).precision - 2 + return "{:.{prec}}".format(num, prec=precision) + + +@safe_repr.register(np.ndarray) +def _(array: np.ndarray) -> str: + return f"{type(array).__name__}([{', '.join(map(safe_repr, array))}])" + + +@safe_repr.register(list) +def _(list_: list) -> str: + return f"[{', '.join(map(safe_repr, list_))}]" + + +@safe_repr.register(tuple) +def _(tuple_: tuple) -> str: + return f"({', '.join(map(safe_repr, tuple_))})" diff --git a/test/FIAT/unit/test_macro.py b/test/FIAT/unit/test_macro.py index 1677d7779..8e367bf20 100644 --- a/test/FIAT/unit/test_macro.py +++ b/test/FIAT/unit/test_macro.py @@ -367,6 +367,27 @@ def test_Ck_basis(cell, order, degree, variant): assert numpy.allclose(local_phis, phis[:, ipts]) +def test_C2_double_alfeld(): + # Construct the quintic C2 spline on the double Alfeld split + # See Section 7.5 of Lai & Schumacher + K = ufc_simplex(2) + DCT = AlfeldSplit(AlfeldSplit(K)) + + degree = 5 + + # C3 on major split facets, C2 elsewhere + order = {} + order[1] = dict.fromkeys(DCT.get_interior_facets(1), 2) + order[1].update(dict.fromkeys(range(3, 6), 3)) + + # C4 at minor split barycenters, C3 at major split barycenter + order[0] = dict.fromkeys(DCT.get_interior_facets(0), 4) + order[0][3] = 3 + + P = CkPolynomialSet(DCT, degree, order=order, variant="bubble") + assert P.get_num_members() == 27 + + def test_distance_to_point_l1(cell): A = AlfeldSplit(cell) dim = A.get_spatial_dimension() diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index f551f4d69..d166ba6ef 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -21,6 +21,7 @@ from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron +from FIAT.reference_element import is_ufc, is_hypercube, default_simplex, flatten_reference_cube, Hypercube point = Point() interval = UFCInterval() @@ -32,6 +33,12 @@ triangle_x_interval = TensorProductCell(triangle, interval) quadrilateral_x_interval = TensorProductCell(quadrilateral, interval) +default_interval = default_simplex(1) +default_triangle = default_simplex(2) +default_tetrahedron = default_simplex(3) +default_interval_x_interval = TensorProductCell(default_interval, default_interval) +default_hypercube = Hypercube(2, default_interval_x_interval) + ufc_tetrahedron_21connectivity = [(0, 1, 2), (0, 3, 4), (1, 3, 5), (2, 4, 5)] ufc_hexahedron_21connectivity = [(0, 1, 4, 5), (2, 3, 6, 7), (0, 2, 8, 9), (1, 3, 10, 11), (4, 6, 8, 10), (5, 7, 9, 11)] @@ -396,6 +403,58 @@ def test_distance_to_point_l1(cell, point, expected): assert np.isclose(cell.distance_to_point_l1(point), expected, rtol=1e-3) +@pytest.mark.parametrize(('cell', 'expected'), + [(interval, True), + (triangle, True), + (quadrilateral, True), + (tetrahedron, True), + (interval_x_interval, True), + (triangle_x_interval, True), + (quadrilateral_x_interval, True), + (hexahedron, True), + (default_interval, False), + (default_triangle, False), + (default_tetrahedron, False), + (default_interval_x_interval, False), + (default_hypercube, False),]) +def test_is_ufc(cell, expected): + assert is_ufc(cell) == expected + + +@pytest.mark.parametrize(('cell', 'expected'), + [(interval, True), + (triangle, False), + (quadrilateral, True), + (tetrahedron, False), + (interval_x_interval, True), + (triangle_x_interval, False), + (quadrilateral_x_interval, True), + (hexahedron, True), + (default_interval, True), + (default_triangle, False), + (default_tetrahedron, False), + (default_interval_x_interval, True), + (default_hypercube, True),]) +def test_is_hypercube(cell, expected): + assert is_hypercube(cell) == expected + + +@pytest.mark.parametrize(('cell'), + [interval, + quadrilateral, + interval_x_interval, + triangle_x_interval, + quadrilateral_x_interval, + hexahedron, + default_interval, + default_interval_x_interval, + default_hypercube]) +def test_flatten_maintains_ufc_status(cell): + ufc_status = is_ufc(cell) + flat_cell = flatten_reference_cube(cell) + assert ufc_status == is_ufc(flat_cell) + + if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__)) diff --git a/test/finat/test_create_fiat_element.py b/test/finat/test_create_fiat_element.py index f99053ab9..6530c8243 100644 --- a/test/finat/test_create_fiat_element.py +++ b/test/finat/test_create_fiat_element.py @@ -59,6 +59,11 @@ def tensor_name(request): ids=lambda x: x.cellname(), scope="module") def ufl_A(request, tensor_name): + if request.param == ufl.quadrilateral: + if tensor_name == "DG": + tensor_name = "DQ" + elif tensor_name == "DG L2": + tensor_name = "DQ L2" return finat.ufl.FiniteElement(tensor_name, request.param, 1) diff --git a/test/finat/test_quadrature.py b/test/finat/test_quadrature.py new file mode 100644 index 000000000..4d1acb9f7 --- /dev/null +++ b/test/finat/test_quadrature.py @@ -0,0 +1,19 @@ +import pytest + +from FIAT import ufc_cell +from finat.quadrature import make_quadrature + + +@pytest.mark.parametrize( + "cell_name", + ["interval", "triangle", "interval * interval", "triangle * interval"] +) +def test_quadrature_rules_are_hashable(cell_name): + ref_cell = ufc_cell(cell_name) + quadrature1 = make_quadrature(ref_cell, 3) + quadrature2 = make_quadrature(ref_cell, 3) + + assert quadrature1 is not quadrature2 + assert hash(quadrature1) == hash(quadrature2) + assert repr(quadrature1) == repr(quadrature2) + assert quadrature1 == quadrature2 From 8771342a48a421024335dbbc827dde287a51d438 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 5 May 2025 17:55:59 +0100 Subject: [PATCH 28/35] First stab at a rewrite of connectivity to respect non UFC ordering --- FIAT/reference_element.py | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index edc902f0d..3bbfba637 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -148,20 +148,37 @@ def __init__(self, shape, vertices, topology): # Given the topology, work out for each entity in the cell, # which other entities it contains. self.sub_entities = {} + self.sub_entities_old = {} for dim, entities in topology.items(): self.sub_entities[dim] = {} + self.sub_entities_old[dim] = {} for e, v in entities.items(): vertices = frozenset(v) sub_entities = [] - + sub_entities_old = [] for dim_, entities_ in topology.items(): for e_, vertices_ in entities_.items(): if vertices.issuperset(vertices_): - sub_entities.append((dim_, e_)) - - # Sort for the sake of determinism and by UFC conventions - self.sub_entities[dim][e] = sorted(sub_entities) + sub_entities_old.append((dim_, e_)) + + for e_ in topology[dim].keys(): + # if dim == self.get_spatial_dimension(): + # # if dimension is whole cell, sub entities have a fixed order + # sub_entities.extend([(dim_, e_)]) + if vertices.issuperset(topology[dim][e_]): + # in order to maintain ordering, extract subentities from vertex numbering + top_val_list = list(topology[dim_].values()) + num_verts_dim = len(list(topology[dim].values())[0]) + for i in range(0, num_verts_dim): + indices = list(range(i, len(top_val_list[0]) + i)) + sub_list = tuple(numpy.take(numpy.array(topology[dim][e_]), indices, mode="wrap").tolist()) + for i, val in topology[dim_].items(): + if set(sub_list) == set(val) and (dim_, i) not in sub_entities: + sub_entities.append((dim_, i)) + + self.sub_entities[dim][e] = list(sub_entities) + self.sub_entities_old[dim][e] = list(sub_entities_old) # Build super-entity dictionary by inverting the sub-entity dictionary self.super_entities = {dim: {entity: [] for entity in topology[dim]} for dim in topology} @@ -183,7 +200,7 @@ def __init__(self, shape, vertices, topology): neighbors = children if dim1 < dim0 else parents d01_entities = tuple(e for d, e in neighbors if d == dim1) self.connectivity[(dim0, dim1)].append(d01_entities) - + print(self.connectivity) # Dictionary with derived cells self._split_cache = {} @@ -1271,7 +1288,11 @@ def compare(self, op, other): if isinstance(other, TensorProductCell): return all(op(a, b) for a, b in zip(self.cells, other.cells)) else: - return op(self, other) + if op == operator.gt or op == operator.lt or operator.ne: + return not op(other, self) + if op == operator.ge or op == operator.le: + return not op(other, self) or operator.eq(self, other) + raise ValueError("Unknown operator in cell comparison") def __gt__(self, other): return self.compare(operator.gt, other) From 2373db4df90c6916af6557dd2a710d6609258f01 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 6 May 2025 12:03:20 +0100 Subject: [PATCH 29/35] refactor, simplify --- FIAT/reference_element.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 3bbfba637..0b551181f 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -162,20 +162,17 @@ def __init__(self, shape, vertices, topology): if vertices.issuperset(vertices_): sub_entities_old.append((dim_, e_)) - for e_ in topology[dim].keys(): - # if dim == self.get_spatial_dimension(): - # # if dimension is whole cell, sub entities have a fixed order - # sub_entities.extend([(dim_, e_)]) - if vertices.issuperset(topology[dim][e_]): - # in order to maintain ordering, extract subentities from vertex numbering - top_val_list = list(topology[dim_].values()) - num_verts_dim = len(list(topology[dim].values())[0]) - for i in range(0, num_verts_dim): - indices = list(range(i, len(top_val_list[0]) + i)) - sub_list = tuple(numpy.take(numpy.array(topology[dim][e_]), indices, mode="wrap").tolist()) - for i, val in topology[dim_].items(): - if set(sub_list) == set(val) and (dim_, i) not in sub_entities: - sub_entities.append((dim_, i)) + # in order to maintain ordering, extract subentities from vertex numbering + entities_of_dim_ = list(entities_.values()) + + from itertools import permutations + # generate all possible sub entities + sub_list = permutations(v, len(entities_of_dim_[0])) + for s in sub_list: + # add the sub entities in the same order as in topology + for i, val in entities_.items(): + if set(s) == set(val) and (dim_, i) not in sub_entities: + sub_entities.append((dim_, i)) self.sub_entities[dim][e] = list(sub_entities) self.sub_entities_old[dim][e] = list(sub_entities_old) From 348280097b0cfedd8320138449f33283cac5c2ac Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 7 May 2025 13:48:55 +0100 Subject: [PATCH 30/35] remove print --- FIAT/reference_element.py | 1 - 1 file changed, 1 deletion(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 0b551181f..cffaf299f 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -197,7 +197,6 @@ def __init__(self, shape, vertices, topology): neighbors = children if dim1 < dim0 else parents d01_entities = tuple(e for d, e in neighbors if d == dim1) self.connectivity[(dim0, dim1)].append(d01_entities) - print(self.connectivity) # Dictionary with derived cells self._split_cache = {} From ef6c59c143b834d067c41af7a42578b5ce6c0fdf Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 16 May 2025 16:32:12 +0100 Subject: [PATCH 31/35] modify to use hasse diagram for fuse sub entities --- FIAT/reference_element.py | 83 +++++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 39 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index cffaf299f..d42de4c67 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -133,7 +133,7 @@ class Cell: """Abstract class for a reference cell. Provides accessors for geometry (vertex coordinates) as well as topology (orderings of vertices that make up edges, faces, etc.""" - def __init__(self, shape, vertices, topology): + def __init__(self, shape, vertices, topology, sub_entities=None): """The constructor takes a shape code, the physical vertices expressed as a list of tuples of numbers, and the topology of a cell. @@ -145,38 +145,42 @@ def __init__(self, shape, vertices, topology): self.vertices = vertices self.topology = topology - # Given the topology, work out for each entity in the cell, - # which other entities it contains. - self.sub_entities = {} - self.sub_entities_old = {} - for dim, entities in topology.items(): - self.sub_entities[dim] = {} - self.sub_entities_old[dim] = {} - - for e, v in entities.items(): - vertices = frozenset(v) - sub_entities = [] - sub_entities_old = [] - for dim_, entities_ in topology.items(): - for e_, vertices_ in entities_.items(): - if vertices.issuperset(vertices_): - sub_entities_old.append((dim_, e_)) - - # in order to maintain ordering, extract subentities from vertex numbering - entities_of_dim_ = list(entities_.values()) - - from itertools import permutations - # generate all possible sub entities - sub_list = permutations(v, len(entities_of_dim_[0])) - for s in sub_list: - # add the sub entities in the same order as in topology - for i, val in entities_.items(): - if set(s) == set(val) and (dim_, i) not in sub_entities: - sub_entities.append((dim_, i)) - - self.sub_entities[dim][e] = list(sub_entities) - self.sub_entities_old[dim][e] = list(sub_entities_old) - + if sub_entities: + self.sub_entities = sub_entities + else: + # If sub entity list not provided + # Given the topology, work out for each entity in the cell, + # which other entities it contains. + self.sub_entities = {} + self.sub_entities_old = {} + for dim, entities in topology.items(): + self.sub_entities[dim] = {} + self.sub_entities_old[dim] = {} + + for e, v in entities.items(): + vertices = frozenset(v) + sub_entities = [] + sub_entities_old = [] + for dim_, entities_ in topology.items(): + for e_, vertices_ in entities_.items(): + if vertices.issuperset(vertices_): + sub_entities_old.append((dim_, e_)) + + # in order to maintain ordering, extract subentities from vertex numbering + entities_of_dim_ = list(entities_.values()) + + from itertools import permutations + # generate all possible sub entities + sub_list = permutations(v, len(entities_of_dim_[0])) + for s in sub_list: + # add the sub entities in the same order as in topology + for i, val in entities_.items(): + if set(s) == set(val) and (dim_, i) not in sub_entities: + sub_entities.append((dim_, i)) + + self.sub_entities[dim][e] = list(sub_entities) + self.sub_entities_old[dim][e] = list(sub_entities_old) + self.sub_entities = self.sub_entities_old # Build super-entity dictionary by inverting the sub-entity dictionary self.super_entities = {dim: {entity: [] for entity in topology[dim]} for dim in topology} for dim0 in topology: @@ -199,6 +203,7 @@ def __init__(self, shape, vertices, topology): self.connectivity[(dim0, dim1)].append(d01_entities) # Dictionary with derived cells self._split_cache = {} + print("connectivity", self.connectivity) def __repr__(self): return f"{type(self).__name__}({self.shape!r}, {safe_repr(self.vertices)}, {self.topology!r})" @@ -396,14 +401,14 @@ class SimplicialComplex(Cell): This consists of list of vertex locations and a topology map defining facets. """ - def __init__(self, shape, vertices, topology): + def __init__(self, shape, vertices, topology, sub_ents=None): # Make sure that every facet has the right number of vertices to be # a simplex. for dim in topology: for entity in topology[dim]: assert len(topology[dim][entity]) == dim + 1 - super().__init__(shape, vertices, topology) + super().__init__(shape, vertices, topology, sub_ents) def compute_normal(self, facet_i, cell=None): """Returns the unit normal vector to facet i of codimension 1.""" @@ -1124,7 +1129,7 @@ def compute_normal(self, i): class TensorProductCell(Cell): """A cell that is the product of FIAT cells.""" - def __init__(self, *cells): + def __init__(self, *cells, sub_entities=None): # Vertices vertices = tuple(tuple(chain(*coords)) for coords in product(*[cell.get_vertices() @@ -1147,7 +1152,7 @@ def __init__(self, *cells): topology[dim] = dict(enumerate(topology[dim][key] for key in sorted(topology[dim]))) - super().__init__(TENSORPRODUCT, vertices, topology) + super().__init__(TENSORPRODUCT, vertices, topology, sub_entities) self.cells = tuple(cells) def __repr__(self): @@ -1375,7 +1380,7 @@ def extrinsic_orientation_permutation_map(self): class Hypercube(Cell): """Abstract class for a reference hypercube""" - def __init__(self, dimension, product): + def __init__(self, dimension, product, sub_entities=None): self.dimension = dimension self.shape = hypercube_shapes[dimension] @@ -1383,7 +1388,7 @@ def __init__(self, dimension, product): verts = product.get_vertices() topology = flatten_entities(pt) - super().__init__(self.shape, verts, topology) + super().__init__(self.shape, verts, topology, sub_entities) self.product = product self.unflattening_map = compute_unflattening_map(pt) From df581a74040fb7cd4918fc9310b0dafe87f3946a Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 16 May 2025 16:53:18 +0100 Subject: [PATCH 32/35] remove print --- FIAT/reference_element.py | 1 - 1 file changed, 1 deletion(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index d42de4c67..82b86670e 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -203,7 +203,6 @@ def __init__(self, shape, vertices, topology, sub_entities=None): self.connectivity[(dim0, dim1)].append(d01_entities) # Dictionary with derived cells self._split_cache = {} - print("connectivity", self.connectivity) def __repr__(self): return f"{type(self).__name__}({self.shape!r}, {safe_repr(self.vertices)}, {self.topology!r})" From dde1c056b4410b81e7f2296df5d287855a3be394 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 23 May 2025 09:41:15 +0100 Subject: [PATCH 33/35] add the ablity to pass facet orientation to facet quadrature rule --- FIAT/quadrature.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index 51e4f45e2..a0837e25b 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -196,10 +196,12 @@ class CollapsedQuadratureTetrahedronRule(CollapsedQuadratureSimplexRule): class FacetQuadratureRule(QuadratureRule): """A quadrature rule on a facet mapped from a reference quadrature rule. """ - def __init__(self, ref_el, entity_dim, entity_id, Q_ref): + def __init__(self, ref_el, entity_dim, entity_id, Q_ref, facet_orientation=None): # Construct the facet of interest facet = ref_el.construct_subelement(entity_dim) facet_topology = ref_el.get_topology()[entity_dim][entity_id] + if facet_orientation: + facet_topology = tuple(facet_orientation.permute(list(facet_topology))) facet.vertices = ref_el.get_vertices_of_subcomplex(facet_topology) # Map reference points and weights on the appropriate facet From ab73328e0f7c3363f2128c0144b88c4ac5552247 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 12 Jun 2025 10:19:56 +0100 Subject: [PATCH 34/35] add clearer error --- finat/element_factory.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/finat/element_factory.py b/finat/element_factory.py index 32725f1e0..dd4617d4a 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -116,6 +116,8 @@ def as_fiat_cell(cell): if not isinstance(cell, ufl.AbstractCell): raise ValueError("Expecting a UFL Cell") if isinstance(cell, ufl.TensorProductCell) and any([hasattr(c, "to_fiat") for c in cell._cells]): + if not all([hasattr(c, "to_fiat") for c in cell._cells]): + raise NotImplementedError("FUSE defined cells cannot be tensor producted with FIAT defined cells") return TensorProductCell(*[c.to_fiat() for c in cell._cells]) try: return cell.to_fiat() From e6657618edde7d691498d2299544264af9662fc9 Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Wed, 1 Oct 2025 09:46:13 +0100 Subject: [PATCH 35/35] Move fuse as_cell dependence out of UFL and into final.ufl (#184) * move as_cell to finat --- .github/workflows/test.yml | 3 +++ finat/element_factory.py | 6 +++--- finat/ufl/__init__.py | 2 +- finat/ufl/finiteelement.py | 4 ++-- finat/ufl/finiteelementbase.py | 17 ++++++++++++++++- finat/ufl/mixedelement.py | 3 +-- finat/ufl/tensorproductelement.py | 4 ++-- pyproject.toml | 1 + test/FIAT/regression/scripts/getreferencerepo | 2 +- 9 files changed, 30 insertions(+), 12 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index b15e2d40b..81872dd9e 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -40,6 +40,9 @@ jobs: run: DATA_REPO_GIT="" python -m pytest --cov=FIAT/ test/FIAT - name: Test FInAT run: DATA_REPO_GIT="" python -m pytest --cov=finat/ --cov=gem/ test/finat + - name: Test FInAT with FUSE + run: | + FIREDRAKE_USE_FUSE=1 DATA_REPO_GIT="" python -m pytest test/finat - name: Coveralls if: ${{ github.repository == 'FEniCS/fiat' && github.head_ref == '' && matrix.os == 'ubuntu-latest' && matrix.python-version == '3.11' }} env: diff --git a/finat/element_factory.py b/finat/element_factory.py index 5e4b4e1a3..961192303 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -112,7 +112,7 @@ def as_fiat_cell(cell): :arg cell: the :class:`ufl.Cell` to convert.""" if isinstance(cell, str): - cell = ufl.as_cell(cell) + cell = finat.ufl.as_cell(cell) if not isinstance(cell, ufl.AbstractCell): raise ValueError("Expecting a UFL Cell") if isinstance(cell, ufl.TensorProductCell) and any([hasattr(c, "to_fiat") for c in cell._cells]): @@ -336,8 +336,8 @@ def convert_fuse_element(element, **kwargs): return finat.fiat_elements.FuseElement(element.triple), set() -hexahedron_tpc = ufl.TensorProductCell(ufl.as_cell("interval"), ufl.as_cell("interval"), ufl.as_cell("interval")) -quadrilateral_tpc = ufl.TensorProductCell(ufl.as_cell("interval"), ufl.as_cell("interval")) +hexahedron_tpc = ufl.TensorProductCell(finat.ufl.as_cell("interval"), finat.ufl.as_cell("interval"), finat.ufl.as_cell("interval")) +quadrilateral_tpc = ufl.TensorProductCell(finat.ufl.as_cell("interval"), finat.ufl.as_cell("interval")) _cache = weakref.WeakKeyDictionary() diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py index fd656db6e..1cb7f3dca 100644 --- a/finat/ufl/__init__.py +++ b/finat/ufl/__init__.py @@ -14,7 +14,7 @@ from finat.ufl.brokenelement import BrokenElement # noqa: F401 from finat.ufl.enrichedelement import EnrichedElement, NodalEnrichedElement # noqa: F401 from finat.ufl.finiteelement import FiniteElement # noqa: F401 -from finat.ufl.finiteelementbase import FiniteElementBase # noqa: F401 +from finat.ufl.finiteelementbase import FiniteElementBase, as_cell # noqa: F401 from finat.ufl.hdivcurl import HCurlElement, HDivElement, WithMapping, HDiv, HCurl # noqa: F401 from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement # noqa: F401 from finat.ufl.restrictedelement import RestrictedElement # noqa: F401 diff --git a/finat/ufl/finiteelement.py b/finat/ufl/finiteelement.py index 2b47685bf..6e9aca628 100644 --- a/finat/ufl/finiteelement.py +++ b/finat/ufl/finiteelement.py @@ -11,9 +11,9 @@ # Modified by Massimiliano Leoni, 2016 # Modified by Matthew Scroggs, 2023 -from ufl.cell import TensorProductCell, as_cell +from ufl.cell import TensorProductCell from finat.ufl.elementlist import canonical_element_description, simplices -from finat.ufl.finiteelementbase import FiniteElementBase +from finat.ufl.finiteelementbase import FiniteElementBase, as_cell from ufl.utils.formatting import istr diff --git a/finat/ufl/finiteelementbase.py b/finat/ufl/finiteelementbase.py index 3a84f908e..1b698b056 100644 --- a/finat/ufl/finiteelementbase.py +++ b/finat/ufl/finiteelementbase.py @@ -15,7 +15,7 @@ from hashlib import md5 from ufl import pullback -from ufl.cell import AbstractCell, as_cell +from ufl.cell import AbstractCell, as_cell as as_cell_ufl from ufl.finiteelement import AbstractFiniteElement from ufl.utils.sequences import product @@ -266,3 +266,18 @@ def pullback(self): return pullback.physical_pullback raise ValueError(f"Unsupported mapping: {self.mapping()}") + + +def as_cell(cell: AbstractCell | str | tuple[AbstractCell, ...]) -> AbstractCell: + import os + try: + import fuse + except ModuleNotFoundError: + fuse = None + if isinstance(cell, str) and bool(os.getenv("FIREDRAKE_USE_FUSE", 0)): + if fuse: + return fuse.constructCellComplex(cell) + else: + raise ModuleNotFoundError("FIREDRAKE_USE_FUSE is active but fuse is not installed") + else: + return as_cell_ufl(cell) diff --git a/finat/ufl/mixedelement.py b/finat/ufl/mixedelement.py index e1dd1ccd5..2963d6d1c 100644 --- a/finat/ufl/mixedelement.py +++ b/finat/ufl/mixedelement.py @@ -13,9 +13,8 @@ import numpy as np -from ufl.cell import as_cell from finat.ufl.finiteelement import FiniteElement -from finat.ufl.finiteelementbase import FiniteElementBase +from finat.ufl.finiteelementbase import FiniteElementBase, as_cell from ufl.permutation import compute_indices from ufl.pullback import MixedPullback, SymmetricPullback from ufl.utils.indexflattening import flatten_multiindex, shape_to_strides, unflatten_index diff --git a/finat/ufl/tensorproductelement.py b/finat/ufl/tensorproductelement.py index 14c3d3037..c64e1f702 100644 --- a/finat/ufl/tensorproductelement.py +++ b/finat/ufl/tensorproductelement.py @@ -13,8 +13,8 @@ from itertools import chain -from ufl.cell import TensorProductCell, as_cell -from finat.ufl.finiteelementbase import FiniteElementBase +from ufl.cell import TensorProductCell +from finat.ufl.finiteelementbase import FiniteElementBase, as_cell from ufl.sobolevspace import DirectionalSobolevSpace diff --git a/pyproject.toml b/pyproject.toml index 749a0f64b..7b798de37 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,6 +12,7 @@ dependencies = [ "scipy", "symengine", "sympy", + "fuse-element @ git+https://github.com/firedrakeproject/fuse.git", ] requires-python = ">=3.10" authors = [ diff --git a/test/FIAT/regression/scripts/getreferencerepo b/test/FIAT/regression/scripts/getreferencerepo index a68e1aa0b..0dc1bfaea 100755 --- a/test/FIAT/regression/scripts/getreferencerepo +++ b/test/FIAT/regression/scripts/getreferencerepo @@ -39,7 +39,7 @@ if [ ! -d "$DATA_DIR" ]; then else pushd $DATA_DIR echo "Found existing reference data repository, pulling new data" - git checkout main + git checkout master if [ $? -ne 0 ]; then echo "Failed to checkout main, check state of reference data directory." exit 1