Macro: Support for symbolic tabulation (#73)

* Macro: Support for symbolic tabulation
* Fix DG(0) on a split
* Macro: apply partition of unity for DG tabulation on interior facets
* Macro: implement PowellSabinSplit
* Macro: Lexicographical sort when merging entity_ids
* Fixes for numpy 2.0

FIAT/barycentric_interpolation.py

@@ -36,6 +36,7 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0):
         numpy.multiply(phi, wts[:, None], out=phi)
         numpy.multiply(1.0 / numpy.sum(phi, axis=0), phi, out=phi)
     phi[phi != phi] = 1.0
+    phi = phi.reshape(-1, *pts.shape[:-1])
 
     phi = sp_simplify(phi)
     results = {(0,): phi}
@@ -63,19 +64,17 @@ def __init__(self, ref_el, pts):
         self.points = pts
         self.x = numpy.array(pts, dtype="d").flatten()
         self.cell_node_map = expansions.compute_cell_point_map(ref_el, pts, unique=False)
-        self.dmats = []
-        self.weights = []
-        self.nodes = []
-        for ibfs in self.cell_node_map:
-            nodes = self.x[ibfs]
-            dmat, wts = make_dmat(nodes)
-            self.dmats.append(dmat)
-            self.weights.append(wts)
-            self.nodes.append(nodes)
+        self.dmats = [None for _ in self.cell_node_map]
+        self.weights = [None for _ in self.cell_node_map]
+        self.nodes = [None for _ in self.cell_node_map]
+        for cell, ibfs in self.cell_node_map.items():
+            self.nodes[cell] = self.x[ibfs]
+            self.dmats[cell], self.weights[cell] = make_dmat(self.nodes[cell])
 
         self.degree = max(len(wts) for wts in self.weights)-1
         self.recurrence_order = self.degree + 1
         super(LagrangeLineExpansionSet, self).__init__(ref_el)
+        self.continuity = None if len(self.x) == sum(len(xk) for xk in self.nodes) else "C0"
 
     def get_num_members(self, n):
         return len(self.points)

FIAT/discontinuous_lagrange.py

@@ -217,7 +217,7 @@ class DiscontinuousLagrange(finite_element.CiarletElement):
     def __new__(cls, ref_el, degree, variant="equispaced"):
         if degree == 0:
             splitting, _ = parse_lagrange_variant(variant, discontinuous=True)
-            if splitting is None:
+            if splitting is None and not ref_el.is_macrocell():
                 # FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size()
                 return P0.P0(ref_el)
         return super(DiscontinuousLagrange, cls).__new__(cls)

FIAT/dual_set.py

@@ -199,8 +199,22 @@ def make_entity_closure_ids(ref_el, entity_ids):
     return entity_closure_ids
 
 
+def lexsort_nodes(ref_el, nodes, entity=None, offset=0):
+    """Sort PointEvaluation nodes in lexicographical ordering."""
+    if len(nodes) > 1 and all(isinstance(node, functional.PointEvaluation) for node in nodes):
+        pts = []
+        for node in nodes:
+            pt, = node.get_point_dict()
+            pts.append(pt)
+        bary = ref_el.compute_barycentric_coordinates(pts, entity=entity)
+        order = list(offset + numpy.lexsort(bary.T))
+    else:
+        order = list(range(offset, offset + len(nodes)))
+    return order
+
+
 def merge_entities(nodes, ref_el, entity_ids, entity_permutations):
-    """Collect DOFs from simplicial complex onto facets of parent cell"""
+    """Collect DOFs from simplicial complex onto facets of parent cell."""
     parent_cell = ref_el.get_parent()
     if parent_cell is None:
         return nodes, ref_el, entity_ids, entity_permutations
@@ -215,5 +229,6 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations):
             cur = len(parent_nodes)
             for child_dim, child_entity in parent_to_children[dim][entity]:
                 parent_nodes.extend(nodes[i] for i in entity_ids[child_dim][child_entity])
-            parent_ids[dim][entity] = list(range(cur, len(parent_nodes)))
+            ids = lexsort_nodes(parent_cell, parent_nodes[cur:], entity=(dim, entity), offset=cur)
+            parent_ids[dim][entity] = ids
     return parent_nodes, parent_cell, parent_ids, parent_permutations

FIAT/expansions.py

@@ -242,14 +242,19 @@ def xi_tetrahedron(eta):
     return xi1, xi2, xi3
 
 
-def apply_mapping(A, b, pts):
+def apply_mapping(A, b, pts, transpose=False):
     """Apply an affine mapping to a column-stacked array of points."""
-    if isinstance(pts, numpy.ndarray) and len(pts.shape) == 2:
-        return numpy.dot(A, pts) + b[:, None]
+    if len(pts) == 0:
+        return pts
+    if transpose:
+        Ax = numpy.dot(pts, A.T)
+        for _ in Ax.shape[:-1]:
+            b = b[None, ...]
     else:
-        m1, m2 = A.shape
-        return [sum((A[i, j] * pts[j] for j in range(m2)), b[i])
-                for i in range(m1)]
+        Ax = numpy.dot(A, pts)
+        for _ in Ax.shape[1:]:
+            b = b[..., None]
+    return Ax + b
 
 
 class ExpansionSet(object):
@@ -275,14 +280,14 @@ def __init__(self, ref_el, scale=None, variant=None):
         self.variant = variant
         sd = ref_el.get_spatial_dimension()
         top = ref_el.get_topology()
-        self.base_ref_el = reference_element.default_simplex(sd)
-        base_verts = self.base_ref_el.get_vertices()
+        base_ref_el = reference_element.default_simplex(sd)
+        base_verts = base_ref_el.get_vertices()
 
         self.affine_mappings = [reference_element.make_affine_mapping(
                                 ref_el.get_vertices_of_subcomplex(top[sd][cell]),
                                 base_verts) for cell in top[sd]]
         if scale is None:
-            scale = math.sqrt(1.0 / self.base_ref_el.volume())
+            scale = math.sqrt(1.0 / base_ref_el.volume())
         self.scale = scale
         self.continuity = "C0" if variant == "bubble" else None
         self.recurrence_order = 2
@@ -355,20 +360,43 @@ def distance(alpha, beta):
     def _tabulate(self, n, pts, order=0):
         """A version of tabulate() that also works for a single point."""
         pts = numpy.asarray(pts)
-        cell_point_map = compute_cell_point_map(self.ref_el, pts)
-        phis = [self._tabulate_on_cell(n, pts[ipts], order, cell=k)
-                for k, ipts in enumerate(cell_point_map)]
+        unique = self.continuity is not None and order == 0
+        cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=unique)
+        phis = {cell: self._tabulate_on_cell(n, pts[ipts], order, cell=cell)
+                for cell, ipts in cell_point_map.items()}
 
-        if len(phis) == 1:
+        if not self.ref_el.is_macrocell():
             return phis[0]
 
+        if pts.dtype == object:
+            # If binning is undefined, scale by the characteristic function of each subcell
+            Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique)
+            for cell, phi in phis.items():
+                for alpha in phi:
+                    phi[alpha] *= Xi[cell]
+        elif not unique:
+            # If binning is not unique, divide by the multiplicity of each point
+            mult = numpy.zeros(pts.shape[:-1])
+            for cell, ipts in cell_point_map.items():
+                mult[ipts] += 1
+            for cell, ipts in cell_point_map.items():
+                phi = phis[cell]
+                for alpha in phi:
+                    phi[alpha] /= mult[None, ipts]
+
+        # Insert subcell tabulations into the corresponding submatrices
+        idx = lambda *args: args if args[1] is Ellipsis else numpy.ix_(*args)
         num_phis = self.get_num_members(n)
         cell_node_map = self.get_cell_node_map(n)
         result = {}
-        for alpha in phis[0]:
-            result[alpha] = numpy.zeros((num_phis,) + pts.shape[:-1])
-            for ibfs, ipts, phi in zip(cell_node_map, cell_point_map, phis):
-                result[alpha][numpy.ix_(ibfs, ipts)] = phi[alpha]
+        base_phi = tuple(phis.values())[0]
+        for alpha in base_phi:
+            dtype = base_phi[alpha].dtype
+            result[alpha] = numpy.zeros((num_phis, *pts.shape[:-1]), dtype=dtype)
+            for cell in cell_point_map:
+                ibfs = cell_node_map[cell]
+                ipts = cell_point_map[cell]
+                result[alpha][idx(ibfs, ipts)] += phis[cell][alpha]
         return result
 
     def tabulate_normal_jumps(self, n, ref_pts, facet, order=0):
@@ -388,29 +416,27 @@ def tabulate_normal_jumps(self, n, ref_pts, facet, order=0):
         cell_node_map = self.get_cell_node_map(n)
 
         num_phis = self.get_num_members(n)
-        results = [numpy.zeros((num_phis,) + (sd,) * (r-1) + pts.shape[:-1])
-                   for r in range(order+1)]
-
-        for k, (ibfs, ipts) in enumerate(zip(cell_node_map, cell_point_map)):
-            if len(ipts) > 0:
-                normal = self.ref_el.compute_normal(facet, cell=k)
-                side = numpy.dot(normal, self.ref_el.compute_normal(facet))
-                phi = self._tabulate_on_cell(n, pts[ipts], order, cell=k)
-                v0 = phi[(0,)*sd]
-                for r in range(order+1):
-                    vr = numpy.zeros((sd,)*r + v0.shape, dtype=v0.dtype)
-                    for index in numpy.ndindex(vr.shape[:r]):
-                        vr[index] = phi[tuple(map(index.count, range(sd)))]
-                    if r > 0:
-                        vr = numpy.tensordot(normal, vr, axes=(0, 0))
-                    vr = vr.transpose((-2, *tuple(range(r-1)), -1))
-
-                    shape_indices = tuple(range(sd) for _ in range(r-1))
-                    indices = numpy.ix_(ibfs, *shape_indices, ipts)
-                    if r == 0 and side < 0:
-                        results[r][indices] -= vr
-                    else:
-                        results[r][indices] += vr
+        results = numpy.zeros((order+1, num_phis, *pts.shape[:-1]))
+
+        for cell in cell_point_map:
+            ipts = cell_point_map[cell]
+            ibfs = cell_node_map[cell]
+            normal = self.ref_el.compute_normal(facet, cell=cell)
+            side = numpy.dot(normal, self.ref_el.compute_normal(facet))
+            phi = self._tabulate_on_cell(n, pts[ipts], order, cell=cell)
+            v0 = phi[(0,)*sd]
+            for r in range(order+1):
+                vr = numpy.zeros((sd,)*r + v0.shape, dtype=v0.dtype)
+                for index in numpy.ndindex(vr.shape[:r]):
+                    vr[index] = phi[tuple(map(index.count, range(sd)))]
+                for _ in range(r):
+                    vr = numpy.tensordot(normal, vr, axes=(0, 0))
+
+                indices = numpy.ix_(ibfs, ipts)
+                if r % 2 == 0 and side < 0:
+                    results[r][indices] -= vr
+                else:
+                    results[r][indices] += vr
         return results
 
     def get_dmats(self, degree, cell=0):
@@ -502,7 +528,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         results = {}
         scale = self.get_scale(cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
         for k in range(order+1):
-            v = numpy.zeros((n + 1, len(xs)), "d")
+            v = numpy.zeros((n + 1, len(xs)), xs.dtype)
             if n >= k:
                 v[k:] = jacobi.eval_jacobi_batch(k, k, n-k, xs)
             for p in range(n + 1):
@@ -597,34 +623,94 @@ def polynomial_cell_node_map(ref_el, n, continuity=None):
     return cell_node_map
 
 
+def compute_l1_distance(ref_el, points, entity=None):
+    """Computes the l1 distances from a simplex entity to a set of points.
+
+    :arg ref_el: a SimplicialComplex.
+    :arg points: an iterable of points.
+    :arg entity: a tuple of the entity dimension and id.
+    :returns: a numpy array with the l1 distance from the entity to each point.
+    """
+    if entity is None:
+        entity = (ref_el.get_spatial_dimension(), 0)
+    dim, entity_id = entity
+    ref_verts = numpy.zeros((dim + 1, dim))
+    ref_verts[range(1, dim + 1), range(dim)] = 1.0
+
+    top = ref_el.get_topology()
+    verts = ref_el.get_vertices_of_subcomplex(top[dim][entity_id])
+    A, b = reference_element.make_affine_mapping(verts, ref_verts)
+    A = numpy.vstack((A, -numpy.sum(A, axis=0)))
+    b = numpy.hstack((b, 1-numpy.sum(b, axis=0)))
+    bary = apply_mapping(A, b, points, transpose=True)
+
+    dist = 0.5 * abs(numpy.sum(abs(bary) - bary, axis=-1))
+    return dist
+
+
 def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12):
     """Maps cells on a simplicial complex to points.
 
     :arg ref_el: a SimplicialComplex.
     :arg pts: an iterable of physical points on the complex.
     :kwarg unique: Are we assigning a unique cell to points on facets?
     :kwarg tol: the absolute tolerance.
-    :returns: a numpy array mapping cell id to points located on that cell.
+    :returns: a dict mapping cell id to points located on that cell.
     """
     top = ref_el.get_topology()
     sd = ref_el.get_spatial_dimension()
     if len(top[sd]) == 1:
-        return (Ellipsis,)
+        return {0: Ellipsis}
 
-    cell_point_map = []
-    ref_vertices = reference_element.ufc_simplex(sd).get_vertices()
-    for cell in top[sd]:
-        verts = ref_el.get_vertices_of_subcomplex(top[sd][cell])
-        A, b = reference_element.make_affine_mapping(verts, ref_vertices)
-        A = numpy.vstack((A, -numpy.sum(A, axis=0)))
-        b = numpy.hstack((b, 1-numpy.sum(b, axis=0)))
-        x = numpy.dot(pts, A.T) + b[None, :]
+    pts = numpy.asarray(pts)
+    if pts.dtype == object:
+        return {cell: Ellipsis for cell in sorted(top[sd])}
 
+    cell_point_map = {}
+    for cell in sorted(top[sd]):
         # Bin points based on l1 distance
-        pts_on_cell = abs(numpy.sum(abs(x) - x, axis=1)) < 2*tol
-        if unique:
-            for other in cell_point_map:
-                pts_on_cell[other] = False
-        ipts = numpy.where(pts_on_cell)[0]
-        cell_point_map.append(ipts)
+        pts_on_cell = compute_l1_distance(ref_el, pts, entity=(sd, cell)) < tol
+        if len(pts_on_cell.shape) == 0:
+            # singleton case
+            if pts_on_cell:
+                cell_point_map[cell] = Ellipsis
+                if unique:
+                    break
+        else:
+            if unique:
+                for other in cell_point_map.values():
+                    pts_on_cell[other] = False
+            ipts = numpy.where(pts_on_cell)[0]
+            if len(ipts) > 0:
+                cell_point_map[cell] = ipts
     return cell_point_map
+
+
+def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12):
+    """Computes the partition of unity functions for each subcell.
+
+    :arg ref_el: a SimplicialComplex.
+    :arg pt: a physical point on the complex.
+    :kwarg unique: Are we assigning a unique cell to points on facets?
+    :kwarg tol: the absolute tolerance.
+    :returns: a list of (weighted) characteristic functions for each subcell.
+    """
+    from sympy import Piecewise
+    sd = ref_el.get_spatial_dimension()
+    top = ref_el.get_topology()
+    # assert singleton point
+    pt = pt.reshape((sd,))
+
+    # Compute characteristic function of each cell
+    otherwise = []
+    masks = []
+    for cell in sorted(top[sd]):
+        inside = compute_l1_distance(ref_el, pt, entity=(sd, cell)) < tol
+        masks.append(Piecewise(*otherwise, (1.0, inside), (0.0, True)))
+        if unique:
+            otherwise.append((0.0, inside))
+    # If the point is on a facet, divide the characteristic function by the facet multiplicity
+    if not unique:
+        mult = sum(masks)
+        masks = [m / mult for m in masks]
+    return masks

FIAT/finite_element.py

@@ -148,9 +148,8 @@ def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref
         with warnings.catch_warnings():
             warnings.filterwarnings("error")
             try:
-                LU, piv = scipy.linalg.lu_factor(V)
-                new_coeffs_flat = scipy.linalg.lu_solve((LU, piv), B, trans=1)
-            except scipy.linalg.LinAlgWarning:
+                new_coeffs_flat = scipy.linalg.solve(V, B, transposed=True)
+            except (scipy.linalg.LinAlgWarning, scipy.linalg.LinAlgError):
                 raise numpy.linalg.LinAlgError("Singular Vandermonde matrix")
 
         new_shp = new_coeffs_flat.shape[:1] + shp[1:]
@@ -247,7 +246,7 @@ def entity_support_dofs(elem, entity_dim):
         points = list(map(entity_transform, quad.get_points()))
 
         # Integrate the square of the basis functions on the facet.
-        vals = numpy.double(elem.tabulate(0, points)[(0,) * dim])
+        vals = elem.tabulate(0, points)[(0,) * dim]
         # Ints contains the square of the basis functions
         # integrated over the facet.
         if elem.value_shape():

FIAT/functional.py

@@ -215,15 +215,15 @@ def __init__(self, ref_el, x, alpha):
     def __call__(self, fn):
         """Evaluate the functional on the function fn. Note that this depends
         on sympy being able to differentiate fn."""
-        x = list(self.deriv_dict.keys())[0]
+        x, = self.deriv_dict
 
-        X = sympy.DeferredVector('x')
-        dX = numpy.asarray([X[i] for i in range(len(x))])
+        X = tuple(sympy.Symbol(f"X[{i}]") for i in range(len(x)))
 
-        dvars = tuple(d for d, a in zip(dX, self.alpha)
+        dvars = tuple(d for d, a in zip(X, self.alpha)
                       for count in range(a))
 
-        return sympy.diff(fn(X), *dvars).evalf(subs=dict(zip(dX, x)))
+        df = sympy.lambdify(X, sympy.diff(fn(X), *dvars))
+        return df(*x)
 
 
 class PointNormalDerivative(Functional):

FIAT/hct.py

@@ -31,7 +31,8 @@ def __init__(self, ref_el, degree, reduced=False):
         k = 2 if reduced else 0
         Q = create_quadrature(rline, degree-1+k)
         qpts = Q.get_points()
-        f_at_qpts = eval_jacobi(0, 0, k, 2.0*qpts - 1)
+        x, = qpts.T
+        f_at_qpts = eval_jacobi(0, 0, k, 2.0*x - 1)
         for e in sorted(top[1]):
             cur = len(nodes)
             nodes.append(IntegralMomentOfNormalDerivative(ref_el, e, Q, f_at_qpts))