Add fuse element infrastructure

FIAT/lagrange.py

@@ -40,7 +40,7 @@ def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=Fal
         entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])]
         if sort_entities:
             # sort the entities by support vertex ids
-            support = [top[dim][entity] for dim, entity in entities]
+            support = [sorted(top[dim][entity]) for dim, entity in entities]
             entities = [entity for verts, entity in sorted(zip(support, entities))]
 
         # make nodes by getting points

FIAT/nedelec.py

@@ -22,7 +22,7 @@ def NedelecSpace2D(ref_el, degree):
         raise Exception("NedelecSpace2D requires 2d reference element")
 
     k = degree - 1
-    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,))
+    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal")
 
     dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1)
     dimPk = expansions.polynomial_dimension(ref_el, k)
@@ -32,7 +32,7 @@ def NedelecSpace2D(ref_el, degree):
                                   for i in range(sd))))
     vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices)
 
-    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1)
+    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal")
     PkH = Pkp1.take(list(range(dimPkm1, dimPk)))
 
     Q = create_quadrature(ref_el, 2 * (k + 1))
@@ -43,8 +43,8 @@ def NedelecSpace2D(ref_el, degree):
 
     CrossX = numpy.dot(numpy.array([[0.0, 1.0], [-1.0, 0.0]]), Qpts.T)
     PkHCrossX_at_Qpts = PkH_at_Qpts[:, None, :] * CrossX[None, :, :]
-    PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), Pkp1_at_Qpts.T)
-
+    PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts),
+                                 Pkp1_at_Qpts.T)
     PkHcrossX = polynomial_set.PolynomialSet(ref_el,
                                              k + 1,
                                              k + 1,

FIAT/polynomial_set.py

@@ -176,15 +176,49 @@ def polynomial_set_union_normalized(A, B):
     not contain any of the same members of the set, as we construct a
     span via SVD.
     """
-    new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0)
+    assert A.get_reference_element() == B.get_reference_element()
+    new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B)
+
+    deg = max(A.get_degree(), B.get_degree())
+    em_deg = max(A.get_embedded_degree(), B.get_embedded_degree())
     coeffs = spanning_basis(new_coeffs)
     return PolynomialSet(A.get_reference_element(),
-                         A.get_degree(),
-                         A.get_embedded_degree(),
+                         deg,
+                         em_deg,
                          A.get_expansion_set(),
                          coeffs)
 
 
+def construct_new_coeffs(ref_el, A, B):
+    """
+    Constructs new coefficients for the set union of A and B
+    If A and B are discontinuous and do not have the same degree the smaller one
+    is extended to match the larger.
+
+    This does not handle the case that A and B have continuity but not the same degree.
+    """
+
+    if A.get_expansion_set().continuity != B.get_expansion_set().continuity:
+        raise ValueError("Continuity of expansion sets does not match.")
+
+    if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None:
+        higher = A if A.get_embedded_degree() > B.get_embedded_degree() else B
+        lower = B if A.get_embedded_degree() > B.get_embedded_degree() else A
+
+        diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1]
+
+        # pad only the 0th axis with the difference in size
+        padding = [(0, 0) for i in range(len(lower.coeffs.shape) - 1)] + [(0, diff)]
+        embedded_coeffs = numpy.pad(lower.coeffs, padding)
+
+        new_coeffs = numpy.concatenate((embedded_coeffs, higher.coeffs), axis=0)
+    elif A.get_embedded_degree() == B.get_embedded_degree():
+        new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0)
+    else:
+        raise NotImplementedError("Extending of coefficients is not implemented for PolynomialSets with continuity and different degrees")
+    return new_coeffs
+
+
 class ONSymTensorPolynomialSet(PolynomialSet):
     """Constructs an orthonormal basis for symmetric-tensor-valued
     polynomials on a reference element.

FIAT/raviart_thomas.py

@@ -20,7 +20,7 @@ def RTSpace(ref_el, degree):
     sd = ref_el.get_spatial_dimension()
 
     k = degree - 1
-    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,))
+    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal")
 
     dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1)
     dimPk = expansions.polynomial_dimension(ref_el, k)
@@ -30,7 +30,7 @@ def RTSpace(ref_el, degree):
                                   for i in range(sd))))
     vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices)
 
-    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1)
+    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal")
     PkH = Pkp1.take(list(range(dimPkm1, dimPk)))
 
     Q = create_quadrature(ref_el, 2 * (k + 1))
@@ -39,12 +39,12 @@ def RTSpace(ref_el, degree):
     # have to work on this through "tabulate" interface
     # first, tabulate PkH at quadrature points
     PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd]
+
     Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd]
 
     x = Qpts.T
     PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :]
     PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T)
-
     PkHx = polynomial_set.PolynomialSet(ref_el,
                                         k,
                                         k + 1,

FIAT/reference_element.py

@@ -518,7 +518,8 @@ def make_points(self, dim, entity_id, order, variant=None, interior=1):
         facet of dimension dim.  Order indicates how many points to
         include in each direction."""
         if dim == 0:
-            return (self.get_vertices()[entity_id], )
+            return (self.get_vertices()[self.get_topology()[dim][entity_id][0]],)
+            # return (self.get_vertices()[entity_id], )
         elif 0 < dim <= self.get_spatial_dimension():
             entity_verts = \
                 self.get_vertices_of_subcomplex(
@@ -1343,6 +1344,93 @@ def extrinsic_orientation_permutation_map(self):
         return a
 
 
+class Hypercube(Cell):
+    """
+    For reference elements based on TensorProductCells"""
+    def __init__(self, product):
+        pt = product.get_topology()
+
+        verts = product.get_vertices()
+        topology = flatten_entities(pt)
+
+        # TODO this should be generalised ?
+        cube_types = {2: QUADRILATERAL, 3: HEXAHEDRON}
+        super().__init__(cube_types[sum(product.get_dimension())], verts, topology)
+
+        self.product = product
+        self.unflattening_map = compute_unflattening_map(pt)
+
+    def get_dimension(self):
+        """Returns the subelement dimension of the cell.  Same as the
+        spatial dimension."""
+        return self.get_spatial_dimension()
+
+    def get_entity_transform(self, dim, entity_i):
+        """Returns a mapping of point coordinates from the
+        `entity_i`-th subentity of dimension `dim` to the cell.
+
+        :arg dim: entity dimension (integer)
+        :arg entity_i: entity number (integer)
+        """
+        d, e = self.unflattening_map[(dim, entity_i)]
+        return self.product.get_entity_transform(d, e)
+
+    def volume(self):
+        """Computes the volume in the appropriate dimensional measure."""
+        return self.product.volume()
+
+    def compute_reference_normal(self, facet_dim, facet_i):
+        """Returns the unit normal in infinity norm to facet_i."""
+        assert facet_dim == 1
+        d, i = self.unflattening_map[(facet_dim, facet_i)]
+        return self.product.compute_reference_normal(d, i)
+
+    def contains_point(self, point, epsilon=0):
+        """Checks if reference cell contains given point
+        (with numerical tolerance as given by the L1 distance (aka 'manhatten',
+        'taxicab' or rectilinear distance) to the cell.
+
+        Parameters
+        ----------
+        point : numpy.ndarray, list or symbolic expression
+            The coordinates of the point.
+        epsilon : float
+            The tolerance for the check.
+
+        Returns
+        -------
+        bool : True if the point is inside the cell, False otherwise.
+
+        """
+        return self.product.contains_point(point, epsilon=epsilon)
+
+    def distance_to_point_l1(self, point, rescale=False):
+        """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear
+        distance) to a point with 0.0 if the point is inside the cell.
+
+        For more information see the docstring for the UFCSimplex method."""
+        return self.product.distance_to_point_l1(point, rescale=rescale)
+
+    def symmetry_group_size(self, dim):
+        return [1, 2, 8][dim]
+
+    def cell_orientation_reflection_map(self):
+        """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``)."""
+        return self.product.cell_orientation_reflection_map()
+
+    def __gt__(self, other):
+        return self.product > other
+
+    def __lt__(self, other):
+        return self.product < other
+
+    def __ge__(self, other):
+        return self.product >= other
+
+    def __le__(self, other):
+        return self.product <= other
+
+
 class UFCQuadrilateral(Cell):
     r"""This is the reference quadrilateral with vertices
     (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0).
@@ -1734,17 +1822,22 @@ def tuple_sum(tree):
 
 
 def is_hypercube(cell):
+    res = False
     if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)):
-        return True
+        res = True
     elif isinstance(cell, TensorProductCell):
-        return reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells])
+        res = reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells])
     else:
-        return False
+        res = False
+    res2 = len(cell.vertices()) == 2 ** cell.get_dimension()
+    assert res == res2
+    return res
 
 
 def flatten_reference_cube(ref_el):
     """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube"""
-    flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()}
+    from finat import as_fiat_cell
+    flattened_cube = {2: as_fiat_cell("quadrilateral"), 3: as_fiat_cell("hexahedron")}
     if numpy.sum(ref_el.get_dimension()) <= 1:
         # Just return point/interval cell arguments
         return ref_el
@@ -1797,6 +1890,7 @@ def compute_unflattening_map(topology_dict):
 
 
 def max_complex(complexes):
+    breakpoint()
     max_cell = max(complexes)
     if all(max_cell >= b for b in complexes):
         return max_cell

finat/__init__.py

@@ -8,7 +8,8 @@
                             Lagrange, Real, Serendipity,                                     # noqa: F401
                             TrimmedSerendipityCurl, TrimmedSerendipityDiv,                   # noqa: F401
                             TrimmedSerendipityEdge, TrimmedSerendipityFace,                  # noqa: F401
-                            Nedelec, NedelecSecondKind, RaviartThomas, Regge)                # noqa: F401
+                            Nedelec, NedelecSecondKind, RaviartThomas, Regge,                # noqa: F401
+                            FuseElement)                                                     # noqa: F401
 
 from .argyris import Argyris                                       # noqa: F401
 from .aw import ArnoldWinther, ArnoldWintherNC                     # noqa: F401

finat/cube.py

@@ -1,12 +1,11 @@
 from __future__ import absolute_import, division, print_function
 
-from FIAT.reference_element import (UFCHexahedron, UFCQuadrilateral,
-                                    compute_unflattening_map, flatten_entities,
+from FIAT.reference_element import (compute_unflattening_map, flatten_entities,
                                     flatten_permutations)
 from FIAT.tensor_product import FlattenedDimensions as FIAT_FlattenedDimensions
 from gem.utils import cached_property
-
 from finat.finiteelementbase import FiniteElementBase
+from finat.element_factory import as_fiat_cell
 
 
 class FlattenedDimensions(FiniteElementBase):
@@ -23,9 +22,9 @@ def __init__(self, element):
     def cell(self):
         dim = self.product.cell.get_spatial_dimension()
         if dim == 2:
-            return UFCQuadrilateral()
+            return as_fiat_cell("quadrilateral")
         elif dim == 3:
-            return UFCHexahedron()
+            return as_fiat_cell("hexahedron")
         else:
             raise NotImplementedError("Cannot guess cell for spatial dimension %s" % dim)
 

finat/element_factory.py

@@ -27,6 +27,7 @@
 import ufl
 
 from FIAT import ufc_cell
+from FIAT.reference_element import TensorProductCell
 
 __all__ = ("as_fiat_cell", "create_base_element",
            "create_element", "supported_elements")
@@ -110,9 +111,16 @@ def as_fiat_cell(cell):
     """Convert a ufl cell to a FIAT cell.
 
     :arg cell: the :class:`ufl.Cell` to convert."""
+    if isinstance(cell, str):
+        cell = ufl.as_cell(cell)
     if not isinstance(cell, ufl.AbstractCell):
         raise ValueError("Expecting a UFL Cell")
-    return ufc_cell(cell)
+    if isinstance(cell, ufl.TensorProductCell) and any([hasattr(c, "to_fiat") for c in cell._cells]):
+        return TensorProductCell(*[c.to_fiat() for c in cell._cells])
+    try:
+        return cell.to_fiat()
+    except AttributeError:
+        return ufc_cell(cell)
 
 
 @singledispatch
@@ -305,8 +313,17 @@ def convert_restrictedelement(element, **kwargs):
     return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps
 
 
-hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval)
-quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval)
+@convert.register(finat.ufl.FuseElement)
+def convert_fuse_element(element, **kwargs):
+    if element.triple.flat:
+        new_elem = element.triple.unflatten()
+        finat_elem, deps = _create_element(new_elem.to_ufl(), **kwargs)
+        return finat.FlattenedDimensions(finat_elem), deps
+    return finat.fiat_elements.FuseElement(element.triple), set()
+
+
+hexahedron_tpc = ufl.TensorProductCell(ufl.as_cell("interval"), ufl.as_cell("interval"), ufl.as_cell("interval"))
+quadrilateral_tpc = ufl.TensorProductCell(ufl.as_cell("interval"), ufl.as_cell("interval"))
 _cache = weakref.WeakKeyDictionary()
 
 

finat/fiat_elements.py

@@ -485,3 +485,8 @@ def __init__(self, cell, degree, variant=None):
 class NedelecSecondKind(VectorFiatElement):
     def __init__(self, cell, degree, variant=None):
         super().__init__(FIAT.NedelecSecondKind(cell, degree, variant=variant))
+
+
+class FuseElement(FiatElement):
+    def __init__(self, triple):
+        super(FuseElement, self).__init__(triple.to_fiat())

finat/ufl/__init__.py

@@ -19,3 +19,4 @@
 from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement  # noqa: F401
 from finat.ufl.restrictedelement import RestrictedElement  # noqa: F401
 from finat.ufl.tensorproductelement import TensorProductElement  # noqa: F401
+from finat.ufl.fuseelement import FuseElement  # noqa: F401