DualSet: BLAS3 optimization for IntegralMomentOfDerivative

FIAT/__init__.py

@@ -62,6 +62,7 @@
 from FIAT.bubble import Bubble, FacetBubble
 from FIAT.hdiv_trace import HDivTrace
 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen
+from FIAT.wuxu import WuXuH3NC, WuXuRobustH3NC
 from FIAT.walkington import Walkington
 from FIAT.histopolation import Histopolation
 from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite  # noqa: F401
@@ -88,6 +89,8 @@
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
                       "Hermite": CubicHermite,
+                      "Nonconforming Wu-Xu": WuXuH3NC,
+                      "Nonconforming Robust Wu-Xu": WuXuRobustH3NC,
                       "Hsieh-Clough-Tocher": HsiehCloughTocher,
                       "QuadraticPowellSabin6": QuadraticPowellSabin6,
                       "QuadraticPowellSabin12": QuadraticPowellSabin12,

FIAT/argyris.py

@@ -7,8 +7,7 @@
 from FIAT import finite_element, polynomial_set, dual_set
 from FIAT.check_format_variant import check_format_variant, parse_quadrature_scheme
 from FIAT.functional import (PointEvaluation, PointDerivative, PointNormalDerivative,
-                             IntegralMoment,
-                             IntegralMomentOfNormalDerivative)
+                             IntegralMoment, IntegralMomentOfDerivative)
 from FIAT.jacobi import eval_jacobi_batch, eval_jacobi_deriv_batch
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.reference_element import TRIANGLE, ufc_simplex
@@ -37,16 +36,17 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
             # edge dofs
             k = degree - 5
             rline = ufc_simplex(1)
-            Q = parse_quadrature_scheme(rline, interpolant_deg+k-1, quad_scheme)
-            x = 2.0 * Q.get_points() - 1.0
-            phis = eval_jacobi_batch(2, 2, k, x)
-            dphis = eval_jacobi_deriv_batch(2, 2, k, x)
+            Q_ref = parse_quadrature_scheme(rline, interpolant_deg+k-1, quad_scheme)
+            x = rline.compute_barycentric_coordinates(Q_ref.get_points())
+            xref = x[:, [1]] - x[:, [0]]
+            phis = eval_jacobi_batch(2, 2, k, xref)
+            dphis = 2*eval_jacobi_deriv_batch(2, 2, k, xref)
             for e in sorted(top[1]):
-                Q_mapped = FacetQuadratureRule(ref_el, 1, e, Q)
-                scale = 2 / Q_mapped.jacobian_determinant()
+                Q = FacetQuadratureRule(ref_el, 1, e, Q_ref, avg=True)
+                n = ref_el.compute_normal(e)
                 cur = len(nodes)
-                nodes.extend(IntegralMomentOfNormalDerivative(ref_el, e, Q, phi) for phi in phis)
-                nodes.extend(IntegralMoment(ref_el, Q_mapped, dphi * scale) for dphi in dphis[1:])
+                nodes.extend(IntegralMomentOfDerivative(ref_el, Q, phi, n) for phi in phis)
+                nodes.extend(IntegralMoment(ref_el, Q, dphi) for dphi in dphis[1:])
                 entity_ids[1][e].extend(range(cur, len(nodes)))
 
             # interior dofs
@@ -55,11 +55,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
                 cell = ref_el.construct_subelement(sd)
                 Q_ref = parse_quadrature_scheme(cell, interpolant_deg + q, quad_scheme)
                 Pq = polynomial_set.ONPolynomialSet(cell, q, scale=1)
-                Phis = Pq.tabulate(Q_ref.get_points())[(0,) * sd]
+                phis = Pq.tabulate(Q_ref.get_points())[(0,) * sd]
                 for entity in sorted(top[sd]):
-                    Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
-                    scale = 1 / Q.jacobian_determinant()
-                    phis = scale * Phis
+                    Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref, avg=True)
                     cur = len(nodes)
                     nodes.extend(IntegralMoment(ref_el, Q, phi) for phi in phis)
                     entity_ids[sd][entity] = list(range(cur, len(nodes)))

FIAT/brezzi_douglas_marini.py

@@ -34,9 +34,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
             Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
             for f in top[sd - 1]:
                 cur = len(nodes)
-                Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref)
-                Jdet = Q.jacobian_determinant()
-                n = ref_el.compute_scaled_normal(f) / Jdet
+                Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref, avg=True)
+                n = ref_el.compute_scaled_normal(f)
                 phis = n[None, :, None] * Pq_at_qpts[:, None, :]
                 nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi)
                              for phi in phis)

FIAT/crouzeix_raviart.py

@@ -38,20 +38,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
                 facet = ref_el.construct_subelement(dim)
                 if dim == 0:
                     Q_facet = parse_quadrature_scheme(facet, degree + interpolant_deg-1, quad_scheme)
-                    Phis = numpy.ones((1, len(Q_facet.pts)))
+                    phis = numpy.ones((1, len(Q_facet.pts)))
                 else:
                     k = degree - 1 if dim == sd-1 else degree - (1+dim)
                     if k < 0:
                         continue
                     Q_facet = parse_quadrature_scheme(facet, k + interpolant_deg, quad_scheme)
                     poly_set = polynomial_set.ONPolynomialSet(facet, k)
-                    Phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim]
+                    phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim]
 
                 for i in sorted(top[dim]):
                     cur = len(nodes)
-                    Q = FacetQuadratureRule(ref_el, dim, i, Q_facet)
-                    scale = 1 / Q.jacobian_determinant()
-                    phis = scale * Phis
+                    Q = FacetQuadratureRule(ref_el, dim, i, Q_facet, avg=True)
                     nodes.extend(functional.IntegralMoment(ref_el, Q, phi) for phi in phis)
                     entity_ids[dim][i].extend(range(cur, len(nodes)))
         else:

FIAT/dual_set.py

@@ -7,6 +7,8 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 import numpy
+from itertools import chain
+from collections import defaultdict
 
 from FIAT import polynomial_set, functional
 from FIAT.reference_element import compute_unflattening_map
@@ -117,36 +119,37 @@ def to_riesz(self, poly_set):
         riesz_shape = (num_nodes, *tshape, num_exp)
         mat = numpy.zeros(riesz_shape, "d")
 
-        pts = set()
-        dpts = set()
-        Qs_to_ells = dict()
-        for i, ell in enumerate(self.nodes):
-            if len(ell.deriv_dict) > 0:
-                dpts.update(ell.deriv_dict.keys())
-                continue
-            if isinstance(ell, functional.IntegralMoment):
-                Q = ell.Q
-            else:
-                Q = None
-                pts.update(ell.pt_dict.keys())
-            if Q in Qs_to_ells:
+        def map_quadratures_to_points(nodes, deriv=False):
+            Qs_to_ells = defaultdict(list)
+            for i, ell in enumerate(nodes):
+                if deriv and len(ell.deriv_dict) == 0:
+                    continue
+                elif not deriv and len(ell.pt_dict) == 0:
+                    continue
+                if isinstance(ell, (functional.IntegralMoment, functional.IntegralMomentOfDerivative)):
+                    Q = ell.Q
+                else:
+                    Q = None
                 Qs_to_ells[Q].append(i)
-            else:
-                Qs_to_ells[Q] = [i]
-
-        Qs_to_pts = {}
-        if len(pts) > 0:
-            Qs_to_pts[None] = tuple(sorted(pts))
-        for Q in Qs_to_ells:
-            if Q is not None:
-                cur_pts = tuple(map(tuple, Q.pts))
+            pts = set()
+            Qs_to_pts = {}
+            for Q in Qs_to_ells:
+                if Q is None:
+                    if deriv:
+                        cur_pts = chain.from_iterable(nodes[i].deriv_dict.keys() for i in Qs_to_ells[None])
+                    else:
+                        cur_pts = chain.from_iterable(nodes[i].pt_dict.keys() for i in Qs_to_ells[None])
+                    cur_pts = tuple(set(cur_pts))
+                else:
+                    cur_pts = tuple(map(tuple, Q.pts))
                 Qs_to_pts[Q] = cur_pts
                 pts.update(cur_pts)
+            pts = list(sorted(pts))
+            return Qs_to_ells, Qs_to_pts, pts
 
         # Now tabulate the function values
-        pts = list(sorted(pts))
+        Qs_to_ells, Qs_to_pts, pts = map_quadratures_to_points(self.nodes)
         expansion_values = numpy.transpose(es.tabulate(ed, pts))
-
         for Q in Qs_to_ells:
             ells = Qs_to_ells[Q]
             cur_pts = Qs_to_pts[Q]
@@ -171,25 +174,35 @@ def to_riesz(self, poly_set):
         # Tabulate the derivative values that are needed
         max_deriv_order = max(ell.max_deriv_order for ell in self.nodes)
         if max_deriv_order > 0:
-            dpts = list(sorted(dpts))
+            Qs_to_ells, Qs_to_pts, pts = map_quadratures_to_points(self.nodes, deriv=True)
             # It's easiest/most efficient to get derivatives of the
             # expansion set through the polynomial set interface.
             # This is creating a short-lived set to do just this.
             coeffs = numpy.eye(num_exp)
             expansion = polynomial_set.PolynomialSet(self.ref_el, ed, ed, es, coeffs)
-            dexpansion_values = expansion.tabulate(dpts, max_deriv_order)
-
-            ells = [k for k, ell in enumerate(self.nodes) if len(ell.deriv_dict) > 0]
-            wshape = (len(ells), *tshape, len(dpts))
-            dwts = {alpha: numpy.zeros(wshape, "d") for alpha in dexpansion_values if sum(alpha) > 0}
-            for i, k in enumerate(ells):
-                ell = self.nodes[k]
-                for pt, wac_list in ell.deriv_dict.items():
-                    j = dpts.index(pt)
-                    for (w, alpha, c) in wac_list:
-                        dwts[alpha][i][c][j] = w
-            for alpha in dwts:
-                mat[ells] += numpy.dot(dwts[alpha], dexpansion_values[alpha].T)
+            dexpansion_values = expansion.tabulate(pts, max_deriv_order)
+            for Q in Qs_to_ells:
+                ells = Qs_to_ells[Q]
+                cur_pts = Qs_to_pts[Q]
+                indices = list(map(pts.index, cur_pts))
+                wshape = (len(ells), *tshape, len(cur_pts))
+                dwts = {alpha: numpy.zeros(wshape, "d") for alpha in dexpansion_values if sum(alpha) > 0}
+                if Q is None:
+                    for i, k in enumerate(ells):
+                        ell = self.nodes[k]
+                        for pt, wac_list in ell.deriv_dict.items():
+                            j = cur_pts.index(pt)
+                            for (w, alpha, c) in wac_list:
+                                dwts[alpha][i][c][j] = w
+                else:
+                    for i, k in enumerate(ells):
+                        ell = self.nodes[k]
+                        for alpha in ell.weights:
+                            dwts[alpha][i][ell.comp][:] = ell.weights[alpha]
+                for alpha in dwts:
+                    wts = dwts[alpha]
+                    expansion_values = dexpansion_values[alpha].T
+                    mat[ells] += numpy.dot(wts, expansion_values[indices])
         return mat
 
     def get_indices(self, restriction_domain, take_closure=True):

FIAT/functional.py

@@ -13,9 +13,10 @@
 # - type information
 
 from itertools import chain
+from collections import defaultdict
 import numpy
 
-from FIAT import polynomial_set, jacobi, quadrature_schemes
+from FIAT import polynomial_set, jacobi, quadrature, quadrature_schemes
 
 
 class Functional(object):
@@ -242,21 +243,16 @@ def __init__(self, ref_el, edge_no, pt, comp=(), shp=()):
 class PointSecondDerivative(Functional):
     """Represents d/ds1 d/ds2 at a point."""
     def __init__(self, ref_el, s1, s2, pt, comp=(), shp=(), nm=None):
+        S = numpy.outer(s1, s2)
         sd = ref_el.get_spatial_dimension()
-        tau = numpy.zeros((sd*(sd+1)//2,))
-
-        alphas = []
-        cur = 0
-        for i in range(sd):
-            for j in range(i, sd):
-                alpha = [0] * sd
+        tau = defaultdict(float)
+        for index in numpy.ndindex(S.shape):
+            alpha = [0] * sd
+            for i in index:
                 alpha[i] += 1
-                alpha[j] += 1
-                alphas.append(tuple(alpha))
-                tau[cur] = s1[i] * s2[j] + (i != j) * s2[i] * s1[j]
-                cur += 1
+            tau[tuple(alpha)] += S[index]
 
-        dpt_dict = {tuple(pt): [(tau[i], alphas[i], comp) for i in range(len(alphas))]}
+        dpt_dict = {tuple(pt): [(tau[alpha], alpha, comp) for alpha in tau]}
 
         super().__init__(ref_el, shp, {}, dpt_dict, nm or "PointSecondDeriv")
 
@@ -320,82 +316,54 @@ def __call__(self, fn):
 
 
 class IntegralMomentOfDerivative(Functional):
-    """Functional giving directional derivative integrated against some function on a facet."""
+    """Functional giving multi-directional derivative integrated against some function.
 
-    def __init__(self, ref_el, s, Q, f_at_qpts, comp=(), shp=()):
-        self.f_at_qpts = f_at_qpts
+    :arg ref_el: a :class:`Cell`.
+    :arg Q: a :class:`QuadratureRule`.
+    :arg f_at_qpts: an array tabulating the function f at the quadrature
+         points.
+    :arg *directions: a list of vectors of directions of differentiation.
+    :arg comp: Optional argument indicating that only a particular
+         component of the input function should be integrated against f
+    :arg shp: Optional argument giving the value shape of input functions.
+    """
+
+    def __init__(self, ref_el, Q, f_at_qpts, *directions, comp=(), shp=(), nm=""):
         self.Q = Q
+        self.f_at_qpts = f_at_qpts
+        self.comp = comp
+
+        S = directions[0]
+        for dj in directions[1:]:
+            S = numpy.outer(S, dj)
 
         sd = ref_el.get_spatial_dimension()
+        tau = defaultdict(float)
+        for index in numpy.ndindex(S.shape):
+            alpha = [0] * sd
+            for i in index:
+                alpha[i] += 1
+            tau[tuple(alpha)] += S[index]
 
         points = Q.get_points()
         weights = numpy.multiply(f_at_qpts, Q.get_weights())
+        self.weights = {alpha: weights*tau[alpha] for alpha in tau}
 
-        alphas = tuple(map(tuple, numpy.eye(sd, dtype=int)))
-        dpt_dict = {tuple(pt): [(wt*s[i], alphas[i], comp) for i in range(sd)]
+        dpt_dict = {tuple(pt): [(wt*tau[alpha], alpha, comp) for alpha in tau]
                     for pt, wt in zip(points, weights)}
 
-        super().__init__(ref_el, shp,
-                         {}, dpt_dict, "IntegralMomentOfDerivative")
+        super().__init__(ref_el, shp, {}, dpt_dict, nm or "IntegralMomentOfDerivative")
 
 
-class IntegralMomentOfNormalDerivative(Functional):
+class IntegralMomentOfNormalDerivative(IntegralMomentOfDerivative):
     """Functional giving normal derivative integrated against some function on a facet."""
 
-    def __init__(self, ref_el, facet_no, Q, f_at_qpts, n=None):
-        if n is None:
-            n = ref_el.compute_normal(facet_no)
-        self.n = n
-        self.f_at_qpts = f_at_qpts
-        self.Q = Q
-
-        sd = ref_el.get_spatial_dimension()
-
+    def __init__(self, ref_el, facet_no, Q_face, f_at_qpts):
+        n = ref_el.compute_normal(facet_no)
         # map points onto facet
-        transform = ref_el.get_entity_transform(sd-1, facet_no)
-        points = transform(Q.get_points())
-        self.dpts = points
-        weights = numpy.multiply(f_at_qpts, Q.get_weights())
-
-        alphas = tuple(map(tuple, numpy.eye(sd, dtype=int)))
-        dpt_dict = {tuple(pt): [(wt*n[i], alphas[i], tuple()) for i in range(sd)]
-                    for pt, wt in zip(points, weights)}
-
-        super().__init__(ref_el, tuple(),
-                         {}, dpt_dict, "IntegralMomentOfNormalDerivative")
-
-
-class IntegralMomentOfBidirectionalDerivative(Functional):
-    """Functional giving second derivative integrated against some function on a facet."""
-
-    def __init__(self, ref_el, Q, f_at_qpts, s1, s2):
-        self.f_at_qpts = f_at_qpts
-        self.Q = Q
-
         sd = ref_el.get_spatial_dimension()
-
-        # map points onto facet
-        points = Q.get_points()
-        self.dpts = points
-        weights = numpy.multiply(f_at_qpts, Q.get_weights())
-
-        tau = numpy.zeros((sd*(sd+1)//2,))
-        alphas = []
-        cur = 0
-        for i in range(sd):
-            for j in range(i, sd):
-                alpha = [0] * sd
-                alpha[i] += 1
-                alpha[j] += 1
-                alphas.append(tuple(alpha))
-                tau[cur] = s1[i] * s2[j] + (i != j) * s2[i] * s1[j]
-                cur += 1
-
-        dpt_dict = {tuple(pt): [(wt*tau[i], alphas[i], tuple()) for i in range(len(alphas))]
-                    for pt, wt in zip(points, weights)}
-
-        super().__init__(ref_el, tuple(),
-                         {}, dpt_dict, "IntegralMomentOfBidirectionalDerivative")
+        Q = quadrature.FacetQuadratureRule(ref_el, sd-1, facet_no, Q_face, avg=True)
+        super().__init__(ref_el, Q, f_at_qpts, n, nm="IntegralMomentOfNormalDerivative")
 
 
 class FrobeniusIntegralMoment(IntegralMoment):

FIAT/gopalakrishnan_lederer_schoberl.py

@@ -26,13 +26,11 @@ def __init__(self, ref_el, degree, quad_scheme=None):
 
             for entity in sorted(top[dim]):
                 cur = len(nodes)
-                Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-                Jdet = Q.jacobian_determinant()
+                Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref, avg=True)
                 for f in ref_el.get_connectivity()[(dim, sd-1)][entity]:
                     normal = ref_el.compute_scaled_normal(f)
                     tangents = ref_el.compute_tangents(sd-1, f)
-                    n = normal / Jdet
-                    nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
+                    nodes.extend(BidirectionalMoment(ref_el, t, normal, Q, phi)
                                  for phi in phis for t in tangents)
                 entity_ids[dim][entity].extend(range(cur, len(nodes)))