Dpc refel coord fix

FIAT/__init__.py

@@ -7,13 +7,15 @@
 # Import finite element classes
 from FIAT.finite_element import FiniteElement, CiarletElement  # noqa: F401
 from FIAT.argyris import Argyris
+from FIAT.bernstein import Bernstein
 from FIAT.bell import Bell
 from FIAT.argyris import QuinticArgyris
 from FIAT.brezzi_douglas_marini import BrezziDouglasMarini
 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini
 from FIAT.discontinuous_lagrange import DiscontinuousLagrange
 from FIAT.discontinuous_taylor import DiscontinuousTaylor
 from FIAT.discontinuous_raviart_thomas import DiscontinuousRaviartThomas
+from FIAT.discontinuous_pc import DPC
 from FIAT.hermite import CubicHermite
 from FIAT.lagrange import Lagrange
 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre
@@ -47,12 +49,14 @@
 # List of supported elements and mapping to element classes
 supported_elements = {"Argyris": Argyris,
                       "Bell": Bell,
+                      "Bernstein": Bernstein,
                       "Brezzi-Douglas-Marini": BrezziDouglasMarini,
                       "Brezzi-Douglas-Fortin-Marini": BrezziDouglasFortinMarini,
                       "Bubble": Bubble,
                       "FacetBubble": FacetBubble,
                       "Crouzeix-Raviart": CrouzeixRaviart,
                       "Discontinuous Lagrange": DiscontinuousLagrange,
+                      "DPC": DPC,
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
                       "Hermite": CubicHermite,

FIAT/bernstein.py

@@ -0,0 +1,199 @@
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2018 Miklós Homolya
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+# License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with FIAT.  If not, see <https://www.gnu.org/licenses/>.
+
+import math
+import numpy
+
+from FIAT.finite_element import FiniteElement
+from FIAT.dual_set import DualSet
+from FIAT.polynomial_set import mis
+
+
+class BernsteinDualSet(DualSet):
+    """The dual basis for Bernstein elements."""
+
+    def __init__(self, ref_el, degree):
+        # Initialise data structures
+        topology = ref_el.get_topology()
+        entity_ids = {dim: {entity_i: []
+                            for entity_i in entities}
+                      for dim, entities in topology.items()}
+
+        # Calculate inverse topology
+        inverse_topology = {vertices: (dim, entity_i)
+                            for dim, entities in topology.items()
+                            for entity_i, vertices in entities.items()}
+
+        # Generate triangular barycentric indices
+        dim = ref_el.get_spatial_dimension()
+        kss = mis(dim + 1, degree)
+
+        # Fill data structures
+        nodes = []
+        for i, ks in enumerate(kss):
+            vertices, = numpy.nonzero(ks)
+            entity_dim, entity_i = inverse_topology[tuple(vertices)]
+            entity_ids[entity_dim][entity_i].append(i)
+
+            # Leave nodes unimplemented for now
+            nodes.append(None)
+
+        super(BernsteinDualSet, self).__init__(nodes, ref_el, entity_ids)
+
+
+class Bernstein(FiniteElement):
+    """A finite element with Bernstein polynomials as basis functions."""
+
+    def __init__(self, ref_el, degree):
+        dual = BernsteinDualSet(ref_el, degree)
+        k = 0  # 0-form
+        super(Bernstein, self).__init__(ref_el, dual, degree, k)
+
+    def degree(self):
+        """The degree of the polynomial space."""
+        return self.get_order()
+
+    def value_shape(self):
+        """The value shape of the finite element functions."""
+        return ()
+
+    def tabulate(self, order, points, entity=None):
+        """Return tabulated values of derivatives up to given order of
+        basis functions at given points.
+
+        :arg order: The maximum order of derivative.
+        :arg points: An iterable of points.
+        :arg entity: Optional (dimension, entity number) pair
+                     indicating which topological entity of the
+                     reference element to tabulate on.  If ``None``,
+                     default cell-wise tabulation is performed.
+        """
+        # Transform points to reference cell coordinates
+        ref_el = self.get_reference_element()
+        if entity is None:
+            entity = (ref_el.get_spatial_dimension(), 0)
+
+        entity_dim, entity_id = entity
+        entity_transform = ref_el.get_entity_transform(entity_dim, entity_id)
+        cell_points = list(map(entity_transform, points))
+
+        # Construct Cartesian to Barycentric coordinate mapping
+        vs = numpy.asarray(ref_el.get_vertices())
+        B2R = numpy.vstack([vs.T, numpy.ones(len(vs))])
+        R2B = numpy.linalg.inv(B2R)
+
+        B = numpy.hstack([cell_points,
+                          numpy.ones((len(cell_points), 1))]).dot(R2B.T)
+
+        # Evaluate everything
+        deg = self.degree()
+        dim = ref_el.get_spatial_dimension()
+        raw_result = {(alpha, i): vec
+                      for i, ks in enumerate(mis(dim + 1, deg))
+                      for o in range(order + 1)
+                      for alpha, vec in bernstein_Dx(B, ks, o, R2B).items()}
+
+        # Rearrange result
+        space_dim = self.space_dimension()
+        dtype = numpy.array(list(raw_result.values())).dtype
+        result = {alpha: numpy.zeros((space_dim, len(cell_points)), dtype=dtype)
+                  for o in range(order + 1)
+                  for alpha in mis(dim, o)}
+        for (alpha, i), vec in raw_result.items():
+            result[alpha][i, :] = vec
+        return result
+
+
+def bernstein_db(points, ks, alpha=None):
+    """Evaluates Bernstein polynomials or its derivative at barycentric
+    points.
+
+    :arg points: array of points in barycentric coordinates
+    :arg ks: exponents defining the Bernstein polynomial
+    :arg alpha: derivative tuple
+
+    :returns: array of Bernstein polynomial values at given points.
+    """
+    points = numpy.asarray(points)
+    ks = numpy.array(tuple(ks))
+
+    N, d_1 = points.shape
+    assert d_1 == len(ks)
+
+    if alpha is None:
+        alpha = numpy.zeros(d_1)
+    else:
+        alpha = numpy.array(tuple(alpha))
+        assert d_1 == len(alpha)
+
+    ls = ks - alpha
+    if any(k < 0 for k in ls):
+        return numpy.zeros(len(points))
+    elif all(k == 0 for k in ls):
+        return numpy.ones(len(points))
+    else:
+        # Calculate coefficient
+        coeff = math.factorial(ks.sum())
+        for k in ls:
+            coeff //= math.factorial(k)
+        return coeff * numpy.prod(points**ls, axis=1)
+
+
+def bernstein_Dx(points, ks, order, R2B):
+    """Evaluates Bernstein polynomials or its derivatives according to
+    reference coordinates.
+
+    :arg points: array of points in BARYCENTRIC COORDINATES
+    :arg ks: exponents defining the Bernstein polynomial
+    :arg alpha: derivative order (returns all derivatives of this
+                specified order)
+    :arg R2B: linear mapping from reference to barycentric coordinates
+
+    :returns: dictionary mapping from derivative tuples to arrays of
+              Bernstein polynomial values at given points.
+    """
+    points = numpy.asarray(points)
+    ks = tuple(ks)
+
+    N, d_1 = points.shape
+    assert d_1 == len(ks)
+
+    # Collect derivatives according to barycentric coordinates
+    Db_map = {alpha: bernstein_db(points, ks, alpha)
+              for alpha in mis(d_1, order)}
+
+    # Arrange derivative tensor (barycentric coordinates)
+    dtype = numpy.array(list(Db_map.values())).dtype
+    Db_shape = (d_1,) * order
+    Db_tensor = numpy.empty(Db_shape + (N,), dtype=dtype)
+    for ds in numpy.ndindex(Db_shape):
+        alpha = [0] * d_1
+        for d in ds:
+            alpha[d] += 1
+        Db_tensor[ds + (slice(None),)] = Db_map[tuple(alpha)]
+
+    # Coordinate transformation: barycentric -> reference
+    result = {}
+    for alpha in mis(d_1 - 1, order):
+        values = Db_tensor
+        for d, k in enumerate(alpha):
+            for _ in range(k):
+                values = R2B[:, d].dot(values)
+        result[alpha] = values
+    return result

FIAT/discontinuous_pc.py

@@ -0,0 +1,119 @@
+# Copyright (C) 2018 Cyrus Cheng (Imperial College London)
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+#
+# Modified by David A. Ham (david.ham@imperial.ac.uk), 2018
+
+from FIAT import finite_element, polynomial_set, dual_set, functional
+from FIAT.P0 import P0Dual
+from FIAT.reference_element import (Point,
+                                    DefaultLine,
+                                    UFCInterval,
+                                    UFCQuadrilateral,
+                                    UFCHexahedron,
+                                    UFCTriangle,
+                                    UFCTetrahedron,
+                                    make_affine_mapping)
+import numpy as np
+
+hypercube_simplex_map = {Point(): Point(),
+                         DefaultLine(): DefaultLine(),
+                         UFCInterval(): UFCInterval(),
+                         UFCQuadrilateral(): UFCTriangle(),
+                         UFCHexahedron(): UFCTetrahedron()}
+
+
+class DPC0(finite_element.CiarletElement):
+    def __init__(self, ref_el):
+        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0)
+        dual = P0Dual(ref_el)
+        degree = 0
+        formdegree = ref_el.get_spatial_dimension()  # n-form
+        super(DPC0, self).__init__(poly_set=poly_set,
+                                   dual=dual,
+                                   order=degree,
+                                   ref_el=ref_el,
+                                   formdegree=formdegree)
+
+
+class DPCDualSet(dual_set.DualSet):
+    """The dual basis for DPC elements.  This class works for
+    hypercubes of any dimension.  Nodes are point evaluation at
+    equispaced points.  This is the discontinuous version where
+    all nodes are topologically associated with the cell itself"""
+
+    def __init__(self, ref_el, degree):
+        entity_ids = {}
+        nodes = []
+
+        # Change coordinates here.
+        # Vertices of the simplex corresponding to the reference element.
+        v_simplex = hypercube_simplex_map[ref_el].get_vertices()
+        # Vertices of the reference element.
+        v_hypercube = ref_el.get_vertices()
+        # For the mapping, first two vertices are unchanged in all dimensions.
+        v_ = [v_hypercube[0], v_hypercube[int(-0.5*len(v_hypercube))]]
+
+        # For dimension 1 upwards,
+        # take the next vertex and map it to the midpoint of the edge/face it belongs to, and shares
+        # with no other points.
+        for d in range(1, ref_el.get_dimension()):
+            v_.append(tuple(np.asarray(v_hypercube[ref_el.get_dimension() - d] +
+                            np.average(np.asarray(v_hypercube[::2]), axis=0))))
+        A, b = make_affine_mapping(v_simplex, tuple(v_))  # Make affine mapping to be used later.
+
+        # make nodes by getting points
+        # need to do this dimension-by-dimension, facet-by-facet
+        top = hypercube_simplex_map[ref_el].get_topology()
+        cube_topology = ref_el.get_topology()
+
+        cur = 0
+        for dim in sorted(top):
+            entity_ids[dim] = {}
+            for entity in sorted(top[dim]):
+                pts_cur = hypercube_simplex_map[ref_el].make_points(dim, entity, degree)
+                pts_cur = [tuple(np.matmul(A, np.array(x)) + b) for x in pts_cur]
+                nodes_cur = [functional.PointEvaluation(ref_el, x)
+                             for x in pts_cur]
+                nnodes_cur = len(nodes_cur)
+                nodes += nodes_cur
+                cur += nnodes_cur
+            for entity in sorted(cube_topology[dim]):
+                entity_ids[dim][entity] = []
+
+        entity_ids[dim][0] = list(range(len(nodes)))
+        super(DPCDualSet, self).__init__(nodes, ref_el, entity_ids)
+
+
+class HigherOrderDPC(finite_element.CiarletElement):
+    """The DPC finite element.  It is what it is."""
+
+    def __init__(self, ref_el, degree):
+        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], degree)
+        dual = DPCDualSet(ref_el, degree)
+        formdegree = ref_el.get_spatial_dimension()  # n-form
+        super(HigherOrderDPC, self).__init__(poly_set=poly_set,
+                                             dual=dual,
+                                             order=degree,
+                                             ref_el=ref_el,
+                                             formdegree=formdegree)
+
+
+def DPC(ref_el, degree):
+    if degree == 0:
+        return DPC0(ref_el)
+    else:
+        return HigherOrderDPC(ref_el, degree)

FIAT/finite_element.py

@@ -120,8 +120,8 @@ class CiarletElement(FiniteElement):
     basis generated from polynomials encoded in a `PolynomialSet`.
     """
 
-    def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine"):
-        ref_el = poly_set.get_reference_element()
+    def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref_el=None):
+        ref_el = ref_el or poly_set.get_reference_element()
         super(CiarletElement, self).__init__(ref_el, dual, order, formdegree, mapping)
 
         # build generalized Vandermonde matrix

FIAT/hdiv_trace.py

@@ -348,19 +348,13 @@ def barycentric_coordinates(points, vertices):
     """
 
     # Form mapping matrix
-    last = np.asarray(vertices[-1])
-    T = np.matrix([np.array(v) - last for v in vertices[:-1]]).T
+    T = (np.asarray(vertices[:-1]) - vertices[-1]).T
     invT = np.linalg.inv(T)
 
-    # Compute the barycentric coordinates for all points
-    coords = []
-    for p in points:
-        y = np.asarray(p) - last
-        bary = invT.dot(y.T)
-        bary = [bary[(0, i)] for i in range(len(y))]
-        bary += [1.0 - sum(bary)]
-        coords.append(bary)
-    return coords
+    points = np.asarray(points)
+    bary = np.einsum("ij,kj->ki", invT, (points - vertices[-1]))
+    last = (1 - bary.sum(axis=1))
+    return np.concatenate([bary, last[..., np.newaxis]], axis=1)
 
 
 def map_from_reference_facet(point, vertices):