From c9f30210d1dc899c07280fe05ca9edd51fbd0583 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 17 Oct 2023 10:12:00 +0100 Subject: [PATCH] refactoring --- FIAT/gauss_legendre.py | 31 +++-- FIAT/gauss_lobatto_legendre.py | 21 +++- FIAT/recursive_points.py | 139 +++++++++++------------ test/unit/test_gauss_legendre.py | 3 +- test/unit/test_gauss_lobatto_legendre.py | 3 +- 5 files changed, 104 insertions(+), 93 deletions(-) diff --git a/FIAT/gauss_legendre.py b/FIAT/gauss_legendre.py index 6b3b4f0fd..7fbff629d 100644 --- a/FIAT/gauss_legendre.py +++ b/FIAT/gauss_legendre.py @@ -8,37 +8,44 @@ # # Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021 -from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.reference_element import POINT, LINE, TRIANGLE, TETRAHEDRON +from FIAT import (finite_element, polynomial_set, dual_set, functional, + quadrature, recursive_points) +from FIAT.reference_element import POINT, LINE, TRIANGLE, TETRAHEDRON, UFCInterval from FIAT.orientation_utils import make_entity_permutations_simplex from FIAT.barycentric_interpolation import LagrangePolynomialSet -from FIAT.recursive_points import make_node_family, recursive_points + + +class GaussLegendrePointSet(recursive_points.RecursivePointSet): + """Recursive point set on simplices based on the Gauss-Legendre points on + the interval""" + def __init__(self): + ref_el = UFCInterval() + lr = quadrature.GaussLegendreQuadratureLineRule + f = lambda n: lr(ref_el, n + 1).pts + super(GaussLegendrePointSet, self).__init__(f) class GaussLegendreDualSet(dual_set.DualSet): - """The dual basis for 1D discontinuous elements with nodes at the + """The dual basis for discontinuous elements with nodes at the (recursive) Gauss-Legendre points.""" - node_family = make_node_family("gl") + point_set = GaussLegendrePointSet() def __init__(self, ref_el, degree): entity_ids = {} entity_permutations = {} - - # make nodes by getting points - # need to do this dimension-by-dimension, facet-by-facet top = ref_el.get_topology() - for dim in sorted(top): entity_ids[dim] = {} entity_permutations[dim] = {} - perms = make_entity_permutations_simplex(dim, degree + 1 if dim == len(top) - 1 else -1) for entity in sorted(top[dim]): entity_ids[dim][entity] = [] - entity_permutations[dim][entity] = perms + entity_permutations[dim][entity] = [] - pts = recursive_points(self.node_family, ref_el.vertices, degree) + # make nodes by getting points + pts = self.point_set.recursive_points(ref_el.get_vertices(), degree) nodes = [functional.PointEvaluation(ref_el, x) for x in pts] entity_ids[dim][0] = list(range(len(nodes))) + entity_permutations[dim][0] = make_entity_permutations_simplex(dim, degree + 1) super(GaussLegendreDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py index 31856da87..294f52585 100644 --- a/FIAT/gauss_lobatto_legendre.py +++ b/FIAT/gauss_lobatto_legendre.py @@ -8,17 +8,26 @@ # # Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021 -from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.reference_element import LINE, TRIANGLE, TETRAHEDRON +from FIAT import (finite_element, polynomial_set, dual_set, functional, + quadrature, recursive_points) +from FIAT.reference_element import LINE, TRIANGLE, TETRAHEDRON, UFCInterval from FIAT.orientation_utils import make_entity_permutations_simplex from FIAT.barycentric_interpolation import LagrangePolynomialSet -from FIAT.recursive_points import make_node_family, make_points + + +class GaussLobattoLegendrePointSet(recursive_points.RecursivePointSet): + + def __init__(self): + ref_el = UFCInterval() + lr = quadrature.GaussLobattoLegendreQuadratureLineRule + f = lambda n: lr(ref_el, n + 1).pts if n else None + super(GaussLobattoLegendrePointSet, self).__init__(f) class GaussLobattoLegendreDualSet(dual_set.DualSet): - """The dual basis for simplex continuous elements with nodes at the + """The dual basis for continuous elements with nodes at the (recursive) Gauss-Lobatto points.""" - node_family = make_node_family("gll") + point_set = GaussLobattoLegendrePointSet() def __init__(self, ref_el, degree): entity_ids = {} @@ -35,7 +44,7 @@ def __init__(self, ref_el, degree): entity_permutations[dim] = {} perms = {0: [0]} if dim == 0 else make_entity_permutations_simplex(dim, degree - dim) for entity in sorted(top[dim]): - pts_cur = make_points(self.node_family, ref_el, dim, entity, degree) + pts_cur = self.point_set.make_points(ref_el, dim, entity, degree) nodes_cur = [functional.PointEvaluation(ref_el, x) for x in pts_cur] nnodes_cur = len(nodes_cur) diff --git a/FIAT/recursive_points.py b/FIAT/recursive_points.py index aa87794c0..3ef47e24f 100644 --- a/FIAT/recursive_points.py +++ b/FIAT/recursive_points.py @@ -6,7 +6,6 @@ # # Written by Pablo D. Brubeck (brubeck@protonmail.com), 2023 -from FIAT import quadrature, reference_element import numpy """ @@ -23,23 +22,22 @@ """ -def multiindex_equal(d, k, interior=0): - """A generator for :math:`d`-tuple multi-indices whose sum is :math:`k`. +def multiindex_equal(d, isum, imin=0): + """A generator for d-tuple multi-indices whose sum is isum and minimum is imin. """ if d <= 0: return - imin = interior - imax = k - (d-1) * imin + imax = isum - (d - 1) * imin if imax < imin: return for i in range(imin, imax): - for a in multiindex_equal(d-1, k-i, interior=imin): + for a in multiindex_equal(d - 1, isum - i, imin=imin): yield a + (i,) - yield (imin,)*(d-1) + (imax,) + yield (imin,) * (d - 1) + (imax,) -class NodeFamily: - """Family of nodes on the unit interval. This class essentially is a +class RecursivePointSet(object): + """Family of points on the unit interval. This class essentially is a lazy-evaluate-and-cache dictionary: the user passes a routine to evaluate entries for unknown keys """ @@ -47,84 +45,83 @@ def __init__(self, f): self._f = f self._cache = {} - def __getitem__(self, key): + def interval_points(self, degree): try: - return self._cache[key] + return self._cache[degree] except KeyError: - x = self._f(key) + x = self._f(degree) if x is None: x_ro = x else: x_ro = numpy.array(x).flatten() x_ro.setflags(write=False) - return self._cache.setdefault(key, x_ro) - + return self._cache.setdefault(degree, x_ro) + + def _recursive(self, alpha): + """The barycentric (d-1)-simplex coordinates for a + multiindex alpha of length d and sum n, based on a 1D node family.""" + d = len(alpha) + n = sum(alpha) + b = numpy.zeros((d,), dtype="d") + xn = self.interval_points(n) + if xn is None: + return b + if d == 2: + b[:] = xn[list(alpha)] + return b + weight = 0.0 + for i in range(d): + w = xn[n - alpha[i]] + alpha_noti = alpha[:i] + alpha[i+1:] + br = self._recursive(alpha_noti) + b[:i] += w * br[:i] + b[i+1:] += w * br[i:] + weight += w + b /= weight + return b -def make_node_family(family): - line = reference_element.UFCInterval() + def recursive_points(self, vertices, order, interior=0): + X = numpy.array(vertices) + get_point = lambda alpha: tuple(numpy.dot(self._recursive(alpha), X)) + return list(map(get_point, multiindex_equal(len(vertices), order, interior))) + + def make_points(self, ref_el, dim, entity_id, order): + """Constructs a lattice of points on the entity_id:th + facet of dimension dim. Order indicates how many points to + include in each direction.""" + if dim == 0: + return (ref_el.get_vertices()[entity_id], ) + elif 0 < dim < ref_el.get_spatial_dimension(): + entity_verts = \ + ref_el.get_vertices_of_subcomplex( + ref_el.get_topology()[dim][entity_id]) + return self.recursive_points(entity_verts, order, 1) + elif dim == ref_el.get_spatial_dimension(): + return self.recursive_points(ref_el.get_vertices(), order, 1) + else: + raise ValueError("illegal dimension") + + +def make_recursive_point_set(family): + from FIAT import quadrature, reference_element + ref_el = reference_element.UFCInterval() if family == "equispaced": f = lambda n: numpy.linspace(0.0, 1.0, n + 1) elif family == "dg_equispaced": f = lambda n: numpy.linspace(1.0/(n+2.0), (n+1.0)/(n+2.0), n + 1) elif family == "gl": lr = quadrature.GaussLegendreQuadratureLineRule - f = lambda n: lr(line, n + 1).pts + f = lambda n: lr(ref_el, n + 1).pts elif family == "gll": lr = quadrature.GaussLobattoLegendreQuadratureLineRule - f = lambda n: lr(line, n + 1).pts if n else None + f = lambda n: lr(ref_el, n + 1).pts if n else None else: raise ValueError("Invalid node family %s" % family) - return NodeFamily(f) - - -def recursive(alpha, family): - """The barycentric d-simplex coordinates for a - multiindex alpha with length n, based on a 1D node family.""" - d = len(alpha) - n = sum(alpha) - b = numpy.zeros((d,), dtype="d") - xn = family[n] - if xn is None: - return b - if d == 2: - b[:] = xn[list(alpha)] - return b - weight = 0.0 - for i in range(d): - w = xn[n - alpha[i]] - alpha_noti = alpha[:i] + alpha[i+1:] - br = recursive(alpha_noti, family) - b[:i] += w * br[:i] - b[i+1:] += w * br[i:] - weight += w - b /= weight - return b - - -def recursive_points(family, vertices, order, interior=0): - X = numpy.array(vertices) - get_point = lambda alpha: tuple(numpy.dot(recursive(alpha, family), X)) - return list(map(get_point, multiindex_equal(len(vertices), order, interior=interior))) - - -def make_points(family, ref_el, dim, entity_id, order): - """Constructs a lattice of points on the entity_id:th - facet of dimension dim. Order indicates how many points to - include in each direction.""" - if dim == 0: - return (ref_el.get_vertices()[entity_id], ) - elif 0 < dim < ref_el.get_spatial_dimension(): - entity_verts = \ - ref_el.get_vertices_of_subcomplex( - ref_el.get_topology()[dim][entity_id]) - return recursive_points(family, entity_verts, order, interior=1) - elif dim == ref_el.get_spatial_dimension(): - return recursive_points(family, ref_el.get_vertices(), order, interior=1) - else: - raise ValueError("illegal dimension") + return RecursivePointSet(f) if __name__ == "__main__": + from FIAT import reference_element from matplotlib import pyplot as plt ref_el = reference_element.ufc_simplex(2) h = numpy.sqrt(3) @@ -140,16 +137,16 @@ def make_points(family, ref_el, dim, entity_id, order): # rule = "equispaced" # dg_rule = "dg_equispaced" - family = make_node_family(rule) - dg_family = make_node_family(dg_rule) + family = make_recursive_point_set(rule) + dg_family = make_recursive_point_set(dg_rule) for d in range(1, 4): - print(make_points(family, reference_element.ufc_simplex(d), d, 0, d)) + print(family.make_points(reference_element.ufc_simplex(d), d, 0, d)) topology = ref_el.get_topology() for dim in topology: for entity in topology[dim]: - pts = make_points(family, ref_el, dim, entity, order) + pts = family.make_points(ref_el, dim, entity, order) if len(pts): x = numpy.array(pts) for r in range(1, 3): @@ -171,7 +168,7 @@ def make_points(family, ref_el, dim, entity_id, order): x0 = sum(x[:d])/d plt.plot(x[:, 0], x[:, 1], "k") - pts = recursive_points(dg_family, ref_el.vertices, order) + pts = dg_family.recursive_points(ref_el.vertices, order) x = numpy.array(pts) for r in range(1, 3): th = r * (2*numpy.pi)/3 diff --git a/test/unit/test_gauss_legendre.py b/test/unit/test_gauss_legendre.py index 62b8cd629..b7ae95d09 100644 --- a/test/unit/test_gauss_legendre.py +++ b/test/unit/test_gauss_legendre.py @@ -45,8 +45,7 @@ def test_gl_basis_values(dim, degree): v = lambda x: sum(x)**test_degree coefs = [n(v) for n in fe.dual.nodes] integral = np.dot(coefs, np.dot(tab, q.wts)) - reference = np.dot([sum(x)**test_degree - for x in q.pts], q.wts) + reference = np.dot([v(x) for x in q.pts], q.wts) assert np.allclose(integral, reference, rtol=1e-14) diff --git a/test/unit/test_gauss_lobatto_legendre.py b/test/unit/test_gauss_lobatto_legendre.py index 7b60631b4..b72a692c5 100644 --- a/test/unit/test_gauss_lobatto_legendre.py +++ b/test/unit/test_gauss_lobatto_legendre.py @@ -45,8 +45,7 @@ def test_gll_basis_values(dim, degree): v = lambda x: sum(x)**test_degree coefs = [n(v) for n in fe.dual.nodes] integral = np.dot(coefs, np.dot(tab, q.wts)) - reference = np.dot([sum(x)**test_degree - for x in q.pts], q.wts) + reference = np.dot([v(x) for x in q.pts], q.wts) assert np.allclose(integral, reference, rtol=1e-14)