Merged in miklos1/ffc-elements (pull request #36)

Bring *all* finite elements to FIAT!

Approved-by: Jan Blechta <blechta@karlin.mff.cuni.cz>

FIAT/__init__.py

@@ -33,7 +33,9 @@
 from FIAT.nodal_enriched import NodalEnrichedElement
 from FIAT.discontinuous import DiscontinuousElement
 from FIAT.hdiv_trace import HDivTrace
+from FIAT.mixed import MixedElement                       # noqa: F401
 from FIAT.restricted import RestrictedElement             # noqa: F401
+from FIAT.quadrature_element import QuadratureElement     # noqa: F401
 
 # Important functionality
 from FIAT.quadrature import make_quadrature               # noqa: F401

FIAT/enriched.py

@@ -17,11 +17,14 @@
 
 from __future__ import absolute_import, print_function, division
 
-import numpy as np
-from copy import copy
+from itertools import chain
+
+import numpy
 
 from FIAT.finite_element import FiniteElement
 from FIAT.dual_set import DualSet
+from FIAT.mixed import concatenate_entity_dofs
+
 
 __all__ = ['EnrichedElement']
 
@@ -35,61 +38,44 @@ class EnrichedElement(FiniteElement):
     """
 
     def __init__(self, *elements):
-        assert len(elements) == 2, "EnrichedElement only implemented for two subelements"
-        A, B = elements
-
         # Firstly, check it makes sense to enrich.  Elements must have:
         # - same reference element
         # - same mapping
         # - same value shape
-        if not A.get_reference_element() == B.get_reference_element():
+        if len(set(e.get_reference_element() for e in elements)) > 1:
             raise ValueError("Elements must be defined on the same reference element")
-        if not A.mapping()[0] == B.mapping()[0]:
+        if len(set(m for e in elements for m in e.mapping())) > 1:
             raise ValueError("Elements must have same mapping")
-        if not A.value_shape() == B.value_shape():
+        if len(set(e.value_shape() for e in elements)) > 1:
             raise ValueError("Elements must have the same value shape")
 
         # order is at least max, possibly more, though getting this
         # right isn't important AFAIK
-        order = max(A.get_order(), B.get_order())
+        order = max(e.get_order() for e in elements)
         # form degree is essentially max (not true for Hdiv/Hcurl,
         # but this will raise an error above anyway).
         # E.g. an H^1 function enriched with an L^2 is now just L^2.
-        if A.get_formdegree() is None or B.get_formdegree() is None:
+        if any(e.get_formdegree() is None for e in elements):
             formdegree = None
         else:
-            formdegree = max(A.get_formdegree(), B.get_formdegree())
+            formdegree = max(e.get_formdegree() for e in elements)
 
         # set up reference element and mapping, following checks above
-        ref_el = A.get_reference_element()
-        mapping = A.mapping()[0]
+        ref_el, = set(e.get_reference_element() for e in elements)
+        mapping, = set(m for e in elements for m in e.mapping())
 
         # set up entity_ids - for each geometric entity, just concatenate
         # the entities of the constituent elements
-        Adofs = A.entity_dofs()
-        Bdofs = B.entity_dofs()
-        offset = A.space_dimension()  # number of entities belonging to A
-        entity_ids = {}
-
-        for ent_dim in Adofs:
-            entity_ids[ent_dim] = {}
-            for ent_dim_index in Adofs[ent_dim]:
-                entlist = copy(Adofs[ent_dim][ent_dim_index])
-                entlist += [c + offset for c in Bdofs[ent_dim][ent_dim_index]]
-                entity_ids[ent_dim][ent_dim_index] = entlist
+        entity_ids = concatenate_entity_dofs(elements)
 
         # set up dual basis - just concatenation
-        nodes = A.dual_basis() + B.dual_basis()
+        nodes = list(chain.from_iterable(e.dual_basis() for e in elements))
         dual = DualSet(nodes, ref_el, entity_ids)
 
         super(EnrichedElement, self).__init__(ref_el, dual, order, formdegree, mapping)
 
-        # Set up constituent elements
-        self.A = A
-        self.B = B
-
         # required degree (for quadrature) is definitely max
-        self.polydegree = max(A.degree(), B.degree())
+        self.polydegree = max(e.degree() for e in elements)
 
         # Store subelements
         self._elements = elements
@@ -116,53 +102,35 @@ def tabulate(self, order, points, entity=None):
         """Return tabulated values of derivatives up to given order of
         basis functions at given points."""
 
-        # Again, simply concatenate at the basis-function level
-        # Number of array dimensions depends on whether the space
-        # is scalar- or vector-valued, so treat these separately.
-
-        Asd = self.A.space_dimension()
-        Bsd = self.B.space_dimension()
-        Atab = self.A.tabulate(order, points, entity)
-        Btab = self.B.tabulate(order, points, entity)
-        npoints = len(points)
-        vs = self.A.value_shape()
-        rank = len(vs)  # scalar: 0, vector: 1
-
-        result = {}
-        for index in Atab:
-            if rank == 0:
-                # scalar valued
-                # Atab[index] and Btab[index] look like
-                # array[basis_fn][point]
-                # We build a new array, which will be the concatenation
-                # of the two subarrays, in the first index.
-
-                temp = np.zeros((Asd + Bsd, npoints),
-                                dtype=Atab[index].dtype)
-                temp[:Asd, :] = Atab[index][:, :]
-                temp[Asd:, :] = Btab[index][:, :]
-
-                result[index] = temp
-            elif rank == 1:
-                # vector valued
-                # Atab[index] and Btab[index] look like
-                # array[basis_fn][x/y/z][point]
-                # We build a new array, which will be the concatenation
-                # of the two subarrays, in the first index.
-
-                temp = np.zeros((Asd + Bsd, vs[0], npoints),
-                                dtype=Atab[index].dtype)
-                temp[:Asd, :, :] = Atab[index][:, :, :]
-                temp[Asd:, :, :] = Btab[index][:, :, :]
-
-                result[index] = temp
-            else:
-                raise NotImplementedError("must be scalar- or vector-valued")
-        return result
+        num_components = numpy.prod(self.value_shape())
+        table_shape = (self.space_dimension(), num_components, len(points))
+
+        table = {}
+        irange = slice(0)
+        for element in self._elements:
+
+            etable = element.tabulate(order, points, entity)
+            irange = slice(irange.stop, irange.stop + element.space_dimension())
+
+            # Insert element table into table
+            for dtuple in etable.keys():
+
+                if dtuple not in table:
+                    if num_components == 1:
+                        table[dtuple] = numpy.zeros((self.space_dimension(), len(points)),
+                                                    dtype=etable[dtuple].dtype)
+                    else:
+                        table[dtuple] = numpy.zeros(table_shape,
+                                                    dtype=etable[dtuple].dtype)
+
+                table[dtuple][irange][:] = etable[dtuple]
+
+        return table
 
     def value_shape(self):
         """Return the value shape of the finite element functions."""
-        return self.A.value_shape()
+        result, = set(e.value_shape() for e in self._elements)
+        return result
 
     def dmats(self):
         """Return dmats: expansion coefficients for basis function

FIAT/finite_element.py

@@ -56,7 +56,7 @@ def get_dual_set(self):
         return self.dual
 
     def get_order(self):
-        """Return the order of the element (may be different from the degree."""
+        """Return the order of the element (may be different from the degree)."""
         return self.order
 
     def dual_basis(self):

FIAT/mixed.py

@@ -0,0 +1,123 @@
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2005-2010 Anders Logg
+#
+# 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/>.
+
+from __future__ import absolute_import, print_function, division
+from six.moves import map
+
+import numpy
+
+from collections import defaultdict
+from operator import add
+from functools import partial
+
+from FIAT.dual_set import DualSet
+from FIAT.finite_element import FiniteElement
+
+
+class MixedElement(FiniteElement):
+    """A FIAT-like representation of a mixed element.
+
+    :arg elements: An iterable of FIAT elements.
+    :arg ref_el: The reference element (optional).
+
+    This object offers tabulation of the concatenated basis function
+    tables along with an entity_dofs dict."""
+    def __init__(self, elements, ref_el=None):
+        elements = tuple(elements)
+
+        cells = set(e.get_reference_element() for e in elements)
+        if ref_el is not None:
+            cells.add(ref_el)
+        ref_el, = cells
+
+        # These functionals are absolutely wrong, they all map from
+        # functions of the wrong shape, and potentially of different
+        # shapes.  However, they are wrong precisely as FFC hacks
+        # expect them to be. :(
+        nodes = [L for e in elements for L in e.dual_basis()]
+
+        entity_dofs = concatenate_entity_dofs(elements)
+
+        dual = DualSet(nodes, ref_el, entity_dofs)
+        super(MixedElement, self).__init__(ref_el, dual, None, mapping=None)
+        self._elements = elements
+
+    def elements(self):
+        return self._elements
+
+    def num_sub_elements(self):
+        return len(self._elements)
+
+    def value_shape(self):
+        return (sum(numpy.prod(e.value_shape(), dtype=int) for e in self.elements()), )
+
+    def mapping(self):
+        return [m for e in self._elements for m in e.mapping()]
+
+    def get_nodal_basis(self):
+        raise NotImplementedError("get_nodal_basis not implemented")
+
+    def tabulate(self, order, points, entity=None):
+        """Tabulate a mixed element by appropriately splatting
+        together the tabulation of the individual elements.
+        """
+        shape = (self.space_dimension(),) + self.value_shape() + (len(points),)
+
+        output = {}
+
+        sub_dims = [0] + list(e.space_dimension() for e in self.elements())
+        sub_cmps = [0] + list(numpy.prod(e.value_shape(), dtype=int)
+                              for e in self.elements())
+        irange = numpy.cumsum(sub_dims)
+        crange = numpy.cumsum(sub_cmps)
+
+        for i, e in enumerate(self.elements()):
+            table = e.tabulate(order, points, entity)
+
+            for d, tab in table.items():
+                try:
+                    arr = output[d]
+                except KeyError:
+                    arr = numpy.zeros(shape, dtype=tab.dtype)
+                    output[d] = arr
+
+                ir = irange[i:i+2]
+                cr = crange[i:i+2]
+                tab = tab.reshape(ir[1] - ir[0], cr[1] - cr[0], -1)
+                arr[slice(*ir), slice(*cr)] = tab
+
+        return output
+
+    def is_nodal(self):
+        """True if primal and dual bases are orthogonal."""
+        return all(e.is_nodal() for e in self._elements)
+
+
+def concatenate_entity_dofs(elements):
+    """Combine the entity_dofs from a list of elements into a combined
+    entity_dof containing the information for the concatenated DoFs of
+    all the elements."""
+    entity_dofs = defaultdict(partial(defaultdict, list))
+    offsets = numpy.cumsum([0] + list(e.space_dimension()
+                                      for e in elements), dtype=int)
+    for i, d in enumerate(e.entity_dofs() for e in elements):
+        for dim, dofs in d.items():
+            for ent, off in dofs.items():
+                entity_dofs[dim][ent] += list(map(partial(add, offsets[i]), off))
+    return entity_dofs

FIAT/quadrature_element.py

@@ -0,0 +1,82 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (C) 2007-2016 Kristian B. Oelgaard
+# Copyright (C) 2017 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 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 Garth N. Wells 2006-2009
+
+from __future__ import absolute_import, print_function, division
+from six.moves import range
+
+# Python modules.
+import numpy
+import six
+
+# FIAT modules.
+from FIAT.dual_set import DualSet
+from FIAT.finite_element import FiniteElement
+from FIAT.functional import PointEvaluation
+
+
+class QuadratureElement(FiniteElement):
+    """A set of quadrature points pretending to be a finite element."""
+
+    def __init__(self, ref_el, points):
+        # Create entity dofs.
+        entity_dofs = {dim: {entity: [] for entity in entities}
+                       for dim, entities in six.iteritems(ref_el.get_topology())}
+        entity_dofs[ref_el.get_dimension()] = {0: list(range(len(points)))}
+
+        # The dual nodes are PointEvaluations at the quadrature points.
+        # FIXME: KBO: Check if this gives expected results for code like evaluate_dof.
+        nodes = [PointEvaluation(ref_el, tuple(point)) for point in points]
+
+        # Construct the dual set
+        dual = DualSet(nodes, ref_el, entity_dofs)
+
+        super(QuadratureElement, self).__init__(ref_el, dual, order=None)
+        self._points = points  # save the quadrature points
+
+    def value_shape(self):
+        "The QuadratureElement is scalar valued"
+        return ()
+
+    def tabulate(self, order, points, entity=None):
+        """Return the identity matrix of size (num_quad_points, num_quad_points),
+        in a format that monomialintegration and monomialtabulation understands."""
+
+        if entity is not None and entity != (self.ref_el.get_dimension(), 0):
+            raise ValueError('QuadratureElement does not "tabulate" on subentities.')
+
+        # Derivatives are not defined on a QuadratureElement
+        if order:
+            raise ValueError("Derivatives are not defined on a QuadratureElement.")
+
+        # Check that incoming points are equal to the quadrature points.
+        if len(points) != len(self._points) or abs(numpy.array(points) - self._points).max() > 1e-12:
+            raise AssertionError("Mismatch of quadrature points!")
+
+        # Return the identity matrix of size len(self._points).
+        values = numpy.eye(len(self._points))
+        dim = self.ref_el.get_spatial_dimension()
+        return {(0,) * dim: values}
+
+    @staticmethod
+    def is_nodal():
+        # No polynomial basis, but still nodal.
+        return True