H inf sobolev space (FEniCS#117)

* Add a Sobolev space spanning C infinity elements. * Put R space in HInf * Test Sobolev space additions. * Use Sobolev space as the restrictions check. * Add sobolev_space method to base finite element class (FEniCS#118) * Add nano swp files to .gitignore * MAke finiteelementbase into an abstract base class and add sobolev_space method * _ * restricted element is L2 * s * make fenicsx tests run * branches * _ not - * remove branch from test config * max * main fenics and ffc Co-authored-by: Matthew Scroggs <[email protected]> Co-authored-by: Garth N. Wells <[email protected]>
firedrakeproject · Nov 1, 2022 · 6257d58 · 6257d58
1 parent d88de25
commit 6257d58
Show file tree

Hide file tree

Showing 14 changed files with 122 additions and 77 deletions.
diff --git a/.gitignore b/.gitignore
@@ -8,3 +8,5 @@ test/.pytest_cache
 *.egg-info
 dist/
 release/
+
+.*.swp
diff --git a/test/test_sobolevspace.py b/test/test_sobolevspace.py
@@ -9,14 +9,16 @@
                  FiniteElement, triangle, interval,
                  quadrilateral, HDiv, HCurl)
 from ufl.sobolevspace import SobolevSpace, DirectionalSobolevSpace
-from ufl import H2, H1, HDiv, HCurl, L2
+from ufl import H2, H1, HDiv, HCurl, L2, HInf
+from math import inf
 
 
 # Construct directional Sobolev spaces, with varying smoothness in
 # spatial coordinates
 H0dx0dy = DirectionalSobolevSpace((0, 0))
 H1dx1dy = DirectionalSobolevSpace((1, 1))
 H2dx2dy = DirectionalSobolevSpace((2, 2))
+Hinfdxinfdy = DirectionalSobolevSpace((inf, inf))
 H1dx = DirectionalSobolevSpace((1, 0))
 H1dy = DirectionalSobolevSpace((0, 1))
 H000 = DirectionalSobolevSpace((0, 0, 0))
@@ -44,6 +46,8 @@ def test_directional_space_relations():
     assert H1dx1dy <= HCurl
     assert H2dx2dy <= H1dx1dy
     assert H2dhH1dz < H1
+    assert Hinfdxinfdy <= HInf
+    assert Hinfdxinfdy < H2dx2dy
     assert H1dz > H2dhH1dz
     assert H1dh < L2
     assert H1dz < L2
@@ -63,7 +67,6 @@ def xtest_contains_mixed():
 
 def test_contains_l2():
     l2_elements = [
-        FiniteElement("Real", triangle, 0),
         FiniteElement("DG", triangle, 0),
         FiniteElement("DG", triangle, 1),
         FiniteElement("DG", triangle, 2),
@@ -132,6 +135,22 @@ def test_contains_h2():
         assert h2_element in H0dx0dy
 
 
+def test_contains_hinf():
+    hinf_elements = [
+        FiniteElement("R", triangle, 0)
+    ]
+    for hinf_element in hinf_elements:
+        assert hinf_element in HInf
+        assert hinf_element in H2
+        assert hinf_element in H2dx2dy
+        assert hinf_element in H1
+        assert hinf_element in H1dx1dy
+        assert hinf_element in HDiv
+        assert hinf_element in HCurl
+        assert hinf_element in L2
+        assert hinf_element in H0dx0dy
+
+
 def test_contains_hdiv():
     hdiv_elements = [
         FiniteElement("RT", triangle, 1),

diff --git a/ufl/__init__.py b/ufl/__init__.py
@@ -68,8 +68,11 @@
     - L2
     - H1
     - H2
+    - HInf
     - HDiv
     - HCurl
+    - HEin
+    - HDivDiv
 
 * Elements::
 
@@ -272,7 +275,7 @@
 )
 
 # Sobolev spaces
-from ufl.sobolevspace import L2, H1, H2, HDiv, HCurl
+from ufl.sobolevspace import L2, H1, H2, HDiv, HCurl, HEin, HDivDiv, HInf
 
 # Finite elements classes
 from ufl.finiteelement import FiniteElementBase, FiniteElement, \
@@ -369,7 +372,7 @@
     'UFLException', 'DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL',
     'as_cell', 'AbstractCell', 'Cell', 'TensorProductCell',
     'as_domain', 'AbstractDomain', 'Mesh', 'MeshView', 'TensorProductMesh',
-    'L2', 'H1', 'H2', 'HCurl', 'HDiv',
+    'L2', 'H1', 'H2', 'HCurl', 'HDiv', 'HInf', 'HEin', 'HDivDiv',
     'SpatialCoordinate',
     'CellVolume', 'CellDiameter', 'Circumradius',
     'MinCellEdgeLength', 'MaxCellEdgeLength',

diff --git a/ufl/algorithms/apply_restrictions.py b/ufl/algorithms/apply_restrictions.py
@@ -14,6 +14,7 @@
 from ufl.corealg.map_dag import map_expr_dag
 from ufl.algorithms.map_integrands import map_integrand_dags
 from ufl.measure import integral_type_to_measure_name
+from ufl.sobolevspace import H1
 
 
 class RestrictionPropagator(MultiFunction):
@@ -128,12 +129,7 @@ def reference_value(self, o):
 
     def coefficient(self, o):
         "Allow coefficients to be unrestricted (apply default if so) if the values are fully continuous across the facet."
-        e = o.ufl_element()
-        d = e.degree()
-        f = e.family()
-        # TODO: Move this choice to the element class?
-        continuous_families = ["Lagrange", "Q", "S"]
-        if (f in continuous_families and d > 0) or f == "Real":
+        if o.ufl_element() in H1:
             # If the coefficient _value_ is _fully_ continuous
             return self._default_restricted(o)  # Must still be computed from one of the sides, we just don't care which
         else:

diff --git a/ufl/finiteelement/brokenelement.py b/ufl/finiteelement/brokenelement.py
@@ -14,7 +14,6 @@ class BrokenElement(FiniteElementBase):
     """The discontinuous version of an existing Finite Element space."""
     def __init__(self, element):
         self._element = element
-        self._repr = "BrokenElement(%s)" % repr(element)
 
         family = "BrokenElement"
         cell = element.cell()
@@ -25,15 +24,18 @@ def __init__(self, element):
         FiniteElementBase.__init__(self, family, cell, degree,
                                    quad_scheme, value_shape, reference_value_shape)
 
+    def __repr__(self):
+        return f"BrokenElement({repr(self._element)})"
+
     def mapping(self):
         return self._element.mapping()
 
     def reconstruct(self, **kwargs):
         return BrokenElement(self._element.reconstruct(**kwargs))
 
     def __str__(self):
-        return "BrokenElement(%s)" % str(self._element)
+        return f"BrokenElement({repr(self._element)})"
 
     def shortstr(self):
-        "Format as string for pretty printing."
-        return "BrokenElement(%s)" % str(self._element.shortstr())
+        """Format as string for pretty printing."""
+        return f"BrokenElement({repr(self._element)})"
diff --git a/ufl/finiteelement/elementlist.py b/ufl/finiteelement/elementlist.py
@@ -16,7 +16,7 @@
 from numpy import asarray
 
 from ufl.log import warning, error
-from ufl.sobolevspace import L2, H1, H2, HDiv, HCurl, HEin, HDivDiv
+from ufl.sobolevspace import L2, H1, H2, HDiv, HCurl, HEin, HDivDiv, HInf
 from ufl.utils.formatting import istr
 from ufl.cell import Cell, TensorProductCell
 
@@ -137,7 +137,7 @@ def show_elements():
 register_element("FacetBubble", "FB", 0, H1, "identity", (2, None), simplices)
 register_element("Quadrature", "Quadrature", 0, L2, "identity", (0, None),
                  any_cell)
-register_element("Real", "R", 0, L2, "identity", (0, 0),
+register_element("Real", "R", 0, HInf, "identity", (0, 0),
                  any_cell + ("TensorProductCell",))
 register_element("Undefined", "U", 0, L2, "identity", (0, None), any_cell)
 register_element("Radau", "Rad", 0, L2, "identity", (0, None), ("interval",))
@@ -460,9 +460,6 @@ def canonical_element_description(family, cell, order, form_degree):
             error('Order "%s" invalid for "%s" finite element.' %
                   (istr(order), family))
 
-    # Override sobolev_space for piecewise constants (TODO: necessary?)
-    if order == 0:
-        sobolev_space = L2
     if value_rank == 2:
         # Tensor valued fundamental elements in HEin have this shape
         if gdim is None or tdim is None:

diff --git a/ufl/finiteelement/enrichedelement.py b/ufl/finiteelement/enrichedelement.py
@@ -61,14 +61,11 @@ def __init__(self, *elements):
                                    quad_scheme, value_shape,
                                    reference_value_shape)
 
-        # Cache repr string
-        self._repr = "%s(%s)" % (class_name, ", ".join(repr(e) for e in self._elements))
-
     def mapping(self):
         return self._elements[0].mapping()
 
     def sobolev_space(self):
-        "Return the underlying Sobolev space."
+        """Return the underlying Sobolev space."""
         elements = [e for e in self._elements]
         if all(e.sobolev_space() == elements[0].sobolev_space()
                for e in elements):
@@ -104,6 +101,9 @@ def is_cellwise_constant(self):
         element is spatially constant over each cell."""
         return all(e.is_cellwise_constant() for e in self._elements)
 
+    def __repr__(self):
+        return "EnrichedElement(" + ", ".join(repr(e) for e in self._elements) + ")"
+
     def __str__(self):
         "Format as string for pretty printing."
         return "<%s>" % " + ".join(str(e) for e in self._elements)
@@ -127,6 +127,9 @@ def is_cellwise_constant(self):
         element is spatially constant over each cell."""
         return False
 
+    def __repr__(self):
+        return "NodalEnrichedElement(" + ", ".join(repr(e) for e in self._elements) + ")"
+
     def __str__(self):
         "Format as string for pretty printing."
         return "<Nodal enriched element(%s)>" % ", ".join(str(e) for e in self._elements)

diff --git a/ufl/finiteelement/finiteelement.py b/ufl/finiteelement/finiteelement.py
@@ -24,10 +24,8 @@
 class FiniteElement(FiniteElementBase):
     "The basic finite element class for all simple finite elements."
     # TODO: Move these to base?
-    __slots__ = ("_short_name",
-                 "_sobolev_space",
-                 "_mapping",
-                 "_variant")
+    __slots__ = ("_short_name", "_sobolev_space",
+                 "_mapping", "_variant", "_repr")
 
     def __new__(cls,
                 family,
@@ -188,11 +186,16 @@ def __init__(self,
             repr(self.family()), repr(self.cell()), repr(self.degree()), quad_str, var_str)
         assert '"' not in self._repr
 
+    def __repr__(self):
+        """Format as string for evaluation as Python object."""
+        return self._repr
+
     def mapping(self):
+        """Return the mapping type for this element ."""
         return self._mapping
 
     def sobolev_space(self):
-        "Return the underlying Sobolev space."
+        """Return the underlying Sobolev space."""
         return self._sobolev_space
 
     def variant(self):

diff --git a/ufl/finiteelement/finiteelementbase.py b/ufl/finiteelement/finiteelementbase.py
@@ -15,24 +15,19 @@
 from ufl.utils.dicts import EmptyDict
 from ufl.log import error
 from ufl.cell import AbstractCell, as_cell
+from abc import ABC, abstractmethod
 
 
-class FiniteElementBase(object):
+class FiniteElementBase(ABC):
     "Base class for all finite elements."
-    __slots__ = ("_family",
-                 "_cell",
-                 "_degree",
-                 "_quad_scheme",
-                 "_value_shape",
-                 "_reference_value_shape",
-                 "_repr",
-                 "__weakref__")
+    __slots__ = ("_family", "_cell", "_degree", "_quad_scheme",
+                 "_value_shape", "_reference_value_shape", "__weakref__")
 
     # TODO: Not all these should be in the base class! In particular
     # family, degree, and quad_scheme do not belong here.
     def __init__(self, family, cell, degree, quad_scheme, value_shape,
                  reference_value_shape):
-        "Initialize basic finite element data."
+        """Initialize basic finite element data."""
         if not isinstance(family, str):
             error("Invalid family type.")
         if not (degree is None or isinstance(degree, (int, tuple))):
@@ -54,12 +49,20 @@ def __init__(self, family, cell, degree, quad_scheme, value_shape,
         self._reference_value_shape = reference_value_shape
         self._quad_scheme = quad_scheme
 
+    @abstractmethod
     def __repr__(self):
-        """Format as string for evaluation as Python object.
+        """Format as string for evaluation as Python object."""
+        pass
 
-        NB! Assumes subclass has assigned its repr string to self._repr.
-        """
-        return self._repr
+    @abstractmethod
+    def sobolev_space(self):
+        """Return the underlying Sobolev space."""
+        pass
+
+    @abstractmethod
+    def mapping(self):
+        """Return the mapping type for this element."""
+        pass
 
     def _ufl_hash_data_(self):
         return repr(self)
@@ -101,10 +104,6 @@ def quadrature_scheme(self):
         "Return quadrature scheme of finite element."
         return self._quad_scheme
 
-    def mapping(self):
-        "Not implemented."
-        error("Missing implementation of mapping().")
-
     def cell(self):
         "Return cell of finite element."
         return self._cell
-Original file line number
+Diff line change
@@ Expand Up / @@ -8,3 +8,5 @@ test/.pytest_cache @@
     *.egg-info
     dist/
     release/
+    .*.swp