diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index f87ab28f1..bce14380f 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -14,7 +14,7 @@ class LagrangeDualSet(dual_set.DualSet): simplices of any dimension. Nodes are point evaluation at equispaced points.""" - def __init__(self, ref_el, degree, family="equi"): + def __init__(self, ref_el, degree, family=None): entity_ids = {} nodes = [] entity_permutations = {} diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index 2843afadc..5e0787be9 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -33,7 +33,7 @@ def get_weights(self): return numpy.array(self.wts) def integrate(self, f): - return sum([w * f(x) for (x, w) in zip(self.pts, self.wts)]) + return sum(w * f(x) for x, w in zip(self.pts, self.wts)) class GaussJacobiQuadratureLineRule(QuadratureRule): diff --git a/FIAT/recursive_points.py b/FIAT/recursive_points.py deleted file mode 100644 index 23a740e89..000000000 --- a/FIAT/recursive_points.py +++ /dev/null @@ -1,103 +0,0 @@ -# Copyright (C) 2015 Imperial College London and others. -# -# This file is part of FIAT (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2023 - -import numpy - -""" -@article{isaac2020recursive, - title={Recursive, parameter-free, explicitly defined interpolation nodes for simplices}, - author={Isaac, Tobin}, - journal={SIAM Journal on Scientific Computing}, - volume={42}, - number={6}, - pages={A4046--A4062}, - year={2020}, - publisher={SIAM} -} -""" - - -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 - imax = isum - (d - 1) * imin - if imax < imin: - return - for i in range(imin, imax): - for a in multiindex_equal(d - 1, isum - i, imin=imin): - yield a + (i,) - yield (imin,) * (d - 1) + (imax,) - - -class RecursivePointSet(object): - """Class to construct recursive points on simplices based on a 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 """ - - def __init__(self, family): - self._family = family - self._cache = {} - - def interval_points(self, degree): - try: - return self._cache[degree] - except KeyError: - x = self._family(degree) - if x is not None: - if not isinstance(x, numpy.ndarray): - x = numpy.array(x) - x = x.reshape((-1,)) - x.setflags(write=False) - return self._cache.setdefault(degree, x) - - 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 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") diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index a06e12f57..5684d14ed 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -69,7 +69,7 @@ def lattice_iter(start, finish, depth): yield jj + [ii] -def make_lattice(verts, n, interior=0, family="equi"): +def make_lattice(verts, n, interior=0, family=None): """Constructs a lattice of points on the simplex defined by verts. For example, the 1:st order lattice will be just the vertices. The optional argument interior specifies how many points from @@ -77,6 +77,8 @@ def make_lattice(verts, n, interior=0, family="equi"): and interior = 0, this function will return the vertices and midpoint, but with interior = 1, it will only return the midpoint.""" + if family is None or family == "equispaced": + family = "equi" family = _decode_family(family) D = len(verts) X = numpy.array(verts) @@ -404,7 +406,7 @@ def compute_face_edge_tangents(self, dim, entity_id): edge_ts.append(vert_coords[dest] - vert_coords[source]) return edge_ts - def make_points(self, dim, entity_id, order, family="equi"): + def make_points(self, dim, entity_id, order, family=None): """Constructs a lattice of points on the entity_id:th facet of dimension dim. Order indicates how many points to include in each direction."""