From b290f6c8c0e6e0e6a4685beb713520629a18ff57 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 11 Nov 2023 08:49:36 +0000 Subject: [PATCH 1/6] style --- FIAT/barycentric_interpolation.py | 4 ++-- FIAT/expansions.py | 7 ++----- FIAT/polynomial_set.py | 11 +++++------ 3 files changed, 9 insertions(+), 13 deletions(-) diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py index faae36954..7f5579ff1 100644 --- a/FIAT/barycentric_interpolation.py +++ b/FIAT/barycentric_interpolation.py @@ -85,14 +85,14 @@ def __init__(self, ref_el, pts, shape=tuple()): if shape == tuple(): coeffs = numpy.eye(num_members) else: - coeffs_shape = tuple([num_members] + list(shape) + [num_exp_functions]) + coeffs_shape = (num_members, *shape, num_exp_functions) coeffs = numpy.zeros(coeffs_shape, "d") # use functional's index_iterator function cur_bf = 0 for idx in index_iterator(shape): n = expansions.polynomial_dimension(ref_el, embedded_degree) for exp_bf in range(n): - cur_idx = tuple([cur_bf] + list(idx) + [exp_bf]) + cur_idx = (cur_bf, *idx, exp_bf) coeffs[cur_idx] = 1.0 cur_bf += 1 diff --git a/FIAT/expansions.py b/FIAT/expansions.py index d45724fbe..eba4c6781 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -209,19 +209,16 @@ def get_dmats(self, degree): def _tabulate_jet(self, degree, pts, order=0): from FIAT.polynomial_set import mis - result = {} D = self.ref_el.get_spatial_dimension() lorder = min(2, order) vals = self._tabulate(degree, numpy.transpose(pts), order=lorder) - base_vals = numpy.array(vals[0]) - base_alpha = (0,) * D - result[base_alpha] = base_vals + result = {(0,) * D: numpy.array(vals[0])} for r in range(1, 1+lorder): vr = numpy.transpose(vals[r], tuple(range(1, r+1)) + (0, r+1)) for indices in product(range(D), repeat=r): alpha = tuple(map(indices.count, range(D))) if alpha not in result: - result[alpha] = vr[indices].reshape(base_vals.shape) + result[alpha] = vr[indices] def distance(alpha, beta): return sum(ai != bi for ai, bi in zip(alpha, beta)) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 8a61aaddc..49314877b 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -132,13 +132,12 @@ def __init__(self, ref_el, degree, shape=tuple()): embedded_degree = degree expansion_set = expansions.ExpansionSet(ref_el) + # set up coefficients if shape == tuple(): coeffs = numpy.eye(num_members) else: - # set up coefficients coeffs_shape = (num_members, *shape, num_exp_functions) coeffs = numpy.zeros(coeffs_shape, "d") - # use functional's index_iterator function cur_bf = 0 for idx in index_iterator(shape): @@ -148,8 +147,8 @@ def __init__(self, ref_el, degree, shape=tuple()): coeffs[cur_idx] = 1.0 cur_bf += 1 - PolynomialSet.__init__(self, ref_el, degree, embedded_degree, - expansion_set, coeffs) + super(ONPolynomialSet, self).__init__(ref_el, degree, embedded_degree, + expansion_set, coeffs) def project(f, U, Q): @@ -240,5 +239,5 @@ def __init__(self, ref_el, degree, size=None): coeffs[cur_bf, j, i, exp_bf] = 1.0 cur_bf += 1 - PolynomialSet.__init__(self, ref_el, degree, embedded_degree, - expansion_set, coeffs) + super(ONSymTensorPolynomialSet, self).__init__(ref_el, degree, embedded_degree, + expansion_set, coeffs) From 452eaaa96bae0983a305df76de681417d51af0c9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 11 Nov 2023 09:03:47 +0000 Subject: [PATCH 2/6] ensure we grab the embedding ExpansionSet --- FIAT/lagrange.py | 2 +- FIAT/nodal_enriched.py | 11 +++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index 540894f5e..d46c589e7 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -48,7 +48,7 @@ class Lagrange(finite_element.CiarletElement): def __init__(self, ref_el, degree, variant="equispaced"): dual = LagrangeDualSet(ref_el, degree, variant=variant) - if ref_el.shape == LINE and variant != "equispaced": + if ref_el.shape == LINE: # In 1D we can use the primal basis as the expansion set, # avoiding any round-off coming from a basis transformation points = [] diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py index cd50a2d6b..45928c993 100644 --- a/FIAT/nodal_enriched.py +++ b/FIAT/nodal_enriched.py @@ -35,16 +35,19 @@ def __init__(self, *elements): "of NodalEnrichedElement are nodal") # Extract common data - ref_el = elements[0].get_reference_element() - expansion_set = elements[0].get_nodal_basis().get_expansion_set() degree = min(e.get_nodal_basis().get_degree() for e in elements) embedded_degree = max(e.get_nodal_basis().get_embedded_degree() for e in elements) order = max(e.get_order() for e in elements) - mapping = elements[0].mapping()[0] formdegree = None if any(e.get_formdegree() is None for e in elements) \ else max(e.get_formdegree() for e in elements) - value_shape = elements[0].value_shape() + # LagrangeExpansionSet set has fixed degree, ensure we grab the embedding one + elem = next(e for e in elements + if e.get_nodal_basis().get_embedded_degree() == embedded_degree) + ref_el = elem.get_reference_element() + expansion_set = elem.get_nodal_basis().get_expansion_set() + mapping = elem.mapping()[0] + value_shape = elem.value_shape() # Sanity check assert all(e.get_nodal_basis().get_reference_element() == From 669d179ec676602c1b08f7e5dd2077e78082a1a1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 11 Nov 2023 21:33:46 +0000 Subject: [PATCH 3/6] obtain lagrange coefficients via interpolation --- FIAT/barycentric_interpolation.py | 14 +++++++++----- FIAT/nodal_enriched.py | 16 ++++++++++++---- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py index 7f5579ff1..17c90c426 100644 --- a/FIAT/barycentric_interpolation.py +++ b/FIAT/barycentric_interpolation.py @@ -29,19 +29,23 @@ class LagrangeLineExpansionSet(expansions.LineExpansionSet): https://doi.org/10.1137/S0036144502417715 Eq. (4.2) & (9.4) """ def __init__(self, ref_el, pts): - self.nodes = numpy.array(pts).flatten() - self.dmat, self.weights = make_dmat(self.nodes) + self.points = pts + self.x = numpy.array(pts).flatten() + self.dmat, self.weights = make_dmat(self.x) super(LagrangeLineExpansionSet, self).__init__(ref_el) def get_num_members(self, n): - return len(self.nodes) + return len(self.points) + + def get_points(self): + return self.points def get_dmats(self, degree): return [self.dmat.T] def tabulate(self, n, pts): - assert n == len(self.nodes)-1 - results = numpy.add.outer(-self.nodes, numpy.array(pts).flatten()) + assert n == len(self.points)-1 + results = numpy.add.outer(-self.x, numpy.array(pts).flatten()) with numpy.errstate(divide='ignore', invalid='ignore'): numpy.reciprocal(results, out=results) numpy.multiply(results, self.weights[:, None], out=results) diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py index 45928c993..010831d5c 100644 --- a/FIAT/nodal_enriched.py +++ b/FIAT/nodal_enriched.py @@ -9,6 +9,7 @@ from FIAT.polynomial_set import PolynomialSet from FIAT.dual_set import DualSet from FIAT.finite_element import CiarletElement +from FIAT.barycentric_interpolation import LagrangeLineExpansionSet __all__ = ['NodalEnrichedElement'] @@ -41,7 +42,7 @@ def __init__(self, *elements): order = max(e.get_order() for e in elements) formdegree = None if any(e.get_formdegree() is None for e in elements) \ else max(e.get_formdegree() for e in elements) - # LagrangeExpansionSet set has fixed degree, ensure we grab the embedding one + # LagrangeExpansionSet has fixed degree, ensure we grab the embedding one elem = next(e for e in elements if e.get_nodal_basis().get_embedded_degree() == embedded_degree) ref_el = elem.get_reference_element() @@ -52,14 +53,21 @@ def __init__(self, *elements): # Sanity check assert all(e.get_nodal_basis().get_reference_element() == ref_el for e in elements) - assert all(type(e.get_nodal_basis().get_expansion_set()) == - type(expansion_set) for e in elements) assert all(e_mapping == mapping for e in elements for e_mapping in e.mapping()) assert all(e.value_shape() == value_shape for e in elements) # Merge polynomial sets - coeffs = _merge_coeffs([e.get_coeffs() for e in elements]) + if isinstance(expansion_set, LagrangeLineExpansionSet): + # Obtain coefficients via interpolation + points = expansion_set.get_points() + coeffs = [e.tabulate(0, points)[(0,)] for e in elements] + else: + assert all(type(e.get_nodal_basis().get_expansion_set()) == + type(expansion_set) for e in elements) + coeffs = [e.get_coeffs() for e in elements] + + coeffs = _merge_coeffs(coeffs) poly_set = PolynomialSet(ref_el, degree, embedded_degree, From 2bce25134ed3bebb8846810ae3b4fedd5be2e05b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 13 Nov 2023 19:01:40 +0000 Subject: [PATCH 4/6] dmats with embedded_degree, add more NodalEnriched tests --- FIAT/polynomial_set.py | 2 +- test/unit/test_fiat.py | 18 +++++++++++------- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 49314877b..3a8d376be 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -95,7 +95,7 @@ def get_embedded_degree(self): def get_dmats(self): if len(self.dmats) == 0: - self.dmats = self.expansion_set.get_dmats(self.degree) + self.dmats = self.expansion_set.get_dmats(self.embedded_degree) return self.dmats def get_reference_element(self): diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 57c0ef5f7..58e2ec1b2 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -379,17 +379,21 @@ def test_empty_bubble(): Bubble(S, 3) -def test_nodal_enriched_implementation(): +@pytest.mark.parametrize('elements', [ + (Lagrange(I, 2), Lagrange(I, 1), Bubble(I, 2)), + (GaussLobattoLegendre(I, 3), Lagrange(I, 1), + RestrictedElement(GaussLobattoLegendre(I, 3), restriction_domain="interior")), + (RaviartThomas(T, 2), + RestrictedElement(RaviartThomas(T, 2), restriction_domain='facet'), + RestrictedElement(RaviartThomas(T, 2), restriction_domain='interior')), +]) +def test_nodal_enriched_implementation(elements): """Following element pair should be the same. This might be fragile to dof reordering but works now. """ - e0 = RaviartThomas(T, 2) - - e1 = NodalEnrichedElement( - RestrictedElement(RaviartThomas(T, 2), restriction_domain='facet'), - RestrictedElement(RaviartThomas(T, 2), restriction_domain='interior') - ) + e0 = elements[0] + e1 = NodalEnrichedElement(*elements[1:]) for attr in ["degree", "get_reference_element", From 355718a1db7f7bc64a28a3e07f800c53ed06086b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 13 Nov 2023 21:51:06 +0000 Subject: [PATCH 5/6] NodallyEnriched poly_set needs to be labelled by the max degree --- FIAT/nodal_enriched.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py index 010831d5c..b98eac831 100644 --- a/FIAT/nodal_enriched.py +++ b/FIAT/nodal_enriched.py @@ -36,7 +36,7 @@ def __init__(self, *elements): "of NodalEnrichedElement are nodal") # Extract common data - degree = min(e.get_nodal_basis().get_degree() for e in elements) + degree = max(e.get_nodal_basis().get_degree() for e in elements) embedded_degree = max(e.get_nodal_basis().get_embedded_degree() for e in elements) order = max(e.get_order() for e in elements) From 32d63e1864c70ddb6d575c60c51eac23d2997172 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 15 Nov 2023 09:59:49 +0000 Subject: [PATCH 6/6] use numpy ndindex --- FIAT/expansions.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index eba4c6781..c17cffa0d 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -9,7 +9,6 @@ import numpy import math -from itertools import product from FIAT import reference_element from FIAT import jacobi @@ -215,7 +214,7 @@ def _tabulate_jet(self, degree, pts, order=0): result = {(0,) * D: numpy.array(vals[0])} for r in range(1, 1+lorder): vr = numpy.transpose(vals[r], tuple(range(1, r+1)) + (0, r+1)) - for indices in product(range(D), repeat=r): + for indices in numpy.ndindex(vr.shape[:r]): alpha = tuple(map(indices.count, range(D))) if alpha not in result: result[alpha] = vr[indices] @@ -258,10 +257,10 @@ def tabulate_jet(self, n, pts, order=1): v0 = vals[(0,)*D] data = [v0] for r in range(1, order+1): - v = numpy.zeros((D,)*r + v0.shape, dtype=v0.dtype) - for index in product(range(D), repeat=r): - v[index] = vals[tuple(map(index.count, range(D)))] - data.append(v.transpose((r, r+1) + tuple(range(r)))) + vr = numpy.zeros((D,)*r + v0.shape, dtype=v0.dtype) + for index in numpy.ndindex(vr.shape[:r]): + vr[index] = vals[tuple(map(index.count, range(D)))] + data.append(vr.transpose((r, r+1) + tuple(range(r)))) return data