Add P1-iso-P2 on interval (#232)

* Add P1 iso P2 on interval

* changelog

* is not ==

* piecewisefunction

* mypy

* 1d piecewise functions

* flake

* isort

Author: mscroggs

CHANGELOG_SINCE_LAST_VERSION.md

@@ -6,3 +6,4 @@
 - Allow element creation on non-default references
 - Corrected Bell element
 - Corrected legendre and lobatto variants of Lagrange
+- Add P1-iso-P2 elemen on interval

README.md

@@ -148,6 +148,7 @@ The reference interval has vertices (0,) and (1,). Its sub-entities are numbered
 - Hermite
 - Lagrange (alternative names: P)
 - Morley-Wang-Xu (alternative names: MWX)
+- P1-iso-P2 (alternative names: P2-iso-P1, iso-P2 P1)
 - serendipity (alternative names: S)
 - Taylor (alternative names: discontinuous Taylor)
 - vector Lagrange (alternative names: vP)

symfem/elements/p1_iso_p2.py

@@ -14,7 +14,45 @@
 from ..geometry import SetOfPoints
 from ..piecewise_functions import PiecewiseFunction
 from ..references import Reference
-from ..symbols import x as x_variables
+
+
+class P1IsoP2Interval(CiarletElement):
+    """P1-iso-P2 finite element on an interval."""
+
+    def __init__(self, reference: Reference, order: int):
+        """Create the element.
+
+        Args:
+            reference: The reference element
+            order: The polynomial order
+        """
+        zero = reference.get_point((sympy.Integer(0), ))
+        half = reference.get_point((sympy.Rational(1, 2), ))
+        one = reference.get_point((sympy.Integer(1), ))
+
+        x = reference.get_inverse_map_to_self()[0]
+        poly: typing.List[FunctionInput] = [
+            PiecewiseFunction({(zero, half): 1 - 2 * x, (half, one): 0}, 1),
+            PiecewiseFunction({(zero, half): 2 * x, (half, one): 2 - 2 * x}, 1),
+            PiecewiseFunction({(zero, half): 0, (half, one): 2 * x - 1}, 1),
+        ]
+
+        dofs: ListOfFunctionals = []
+        for v_n, v in enumerate(reference.vertices):
+            dofs.append(PointEvaluation(reference, v, entity=(0, v_n)))
+        entity = reference.sub_entity(1, 0)
+        dofs.append(PointEvaluation(reference, entity.midpoint(), entity=(1, 0)))
+
+        super().__init__(
+            reference, order, poly, dofs, reference.tdim, 1
+        )
+
+    names = ["P1-iso-P2", "P2-iso-P1", "iso-P2 P1"]
+    references = ["interval"]
+    min_order = 1
+    max_order = 1
+    continuity = "C0"
+    last_updated = "2023.08"
 
 
 class P1IsoP2Tri(CiarletElement):
@@ -37,10 +75,10 @@ def __init__(self, reference: Reference, order: int):
             ((zero, half), (half, half), (half, zero)),
         ]
         poly: typing.List[FunctionInput] = []
-        x = x_variables[0]
-        y = x_variables[1]
-        c = 1 - x - y
         invmap = reference.get_inverse_map_to_self()
+        x = invmap[0]
+        y = invmap[1]
+        c = 1 - x - y
         for pieces in [
             {0: 2 * c - 1},
             {1: 2 * x - 1},
@@ -52,7 +90,7 @@ def __init__(self, reference: Reference, order: int):
             poly.append(PiecewiseFunction({
                 tuple(
                     reference.get_point(pt) for pt in q
-                ): pieces[i].subs(x, invmap[0]).subs(y, invmap[1]) if i in pieces else 0
+                ): pieces[i] if i in pieces else 0
                 for i, q in enumerate(tris)
             }, 2))
 
@@ -95,9 +133,9 @@ def __init__(self, reference: Reference, order: int):
             ((half, half), (one, half), (half, one), (one, one)),
         ]
         poly: typing.List[FunctionInput] = []
-        x = x_variables[0]
-        y = x_variables[1]
         invmap = reference.get_inverse_map_to_self()
+        x = invmap[0]
+        y = invmap[1]
         for pieces in [
             {0: (1 - 2 * x) * (1 - 2 * y)},
             {1: (2 * x - 1) * (1 - 2 * y)},
@@ -112,7 +150,7 @@ def __init__(self, reference: Reference, order: int):
             poly.append(PiecewiseFunction({
                 tuple(
                     reference.get_point(pt) for pt in q
-                ): pieces[i].subs(x, invmap[0]).subs(y, invmap[1]) if i in pieces else 0
+                ): pieces[i] if i in pieces else 0
                 for i, q in enumerate(quads)}, 2))
 
         dofs: ListOfFunctionals = []

symfem/finite_element.py

@@ -357,7 +357,10 @@ def get_piece(f, point):
                         return tuple(get_piece(g, point) for g in f)
                     return f
 
-                if self.reference.tdim == 2:
+                if self.reference.tdim == 1:
+                    f = get_piece(f, (0, ))
+                    g = get_piece(g, (0, ))
+                elif self.reference.tdim == 2:
                     f = get_piece(f, (0, sympy.Rational(1, 2)))
                     g = get_piece(g, (0, sympy.Rational(1, 2)))
                 elif self.reference.tdim == 3:

symfem/geometry.py

@@ -87,6 +87,25 @@ def _vdot(v: PointType, w: PointType) -> sympy.core.expr.Expr:
     return out
 
 
+def point_in_interval(point: PointType, interval: SetOfPoints) -> bool:
+    """Check if a point is inside an interval.
+
+    Args:
+        point: The point
+        interval: The vertices of the interval
+
+    Returns:
+        Is the point inside the interval?
+    """
+    v = _vsub(point, interval[0])[0] / _vsub(interval[1], interval[0])[0]
+
+    if isinstance(v, sympy.Float) and _is_close(v, 0):
+        v = sympy.Integer(0)
+    if isinstance(v, sympy.Float) and _is_close(v, 1):
+        v = sympy.Integer(1)
+    return 0 <= v <= 1
+
+
 def point_in_triangle(point: PointType, triangle: SetOfPoints) -> bool:
     """Check if a point is inside a triangle.
 

symfem/piecewise_functions.py

@@ -9,7 +9,8 @@
 from .functions import (AnyFunction, FunctionInput, ScalarFunction, SympyFormat, ValuesToSubstitute,
                         VectorFunction, _to_sympy_format, parse_function_input)
 from .geometry import (PointType, SetOfPoints, SetOfPointsInput, parse_set_of_points_input,
-                       point_in_quadrilateral, point_in_tetrahedron, point_in_triangle)
+                       point_in_interval, point_in_quadrilateral, point_in_tetrahedron,
+                       point_in_triangle)
 from .references import Reference
 from .symbols import AxisVariables, AxisVariablesNotSingle, t, x
 
@@ -113,6 +114,10 @@ def get_piece(self, point: PointType) -> AnyFunction:
         Returns:
             The piece of the function that is valid at that point
         """
+        if self.tdim == 1:
+            for cell, value in self._pieces.items():
+                if point_in_interval(point[:1], cell):
+                    return value
         if self.tdim == 2:
             for cell, value in self._pieces.items():
                 if len(cell) == 3:

symfem/references.py

@@ -220,7 +220,7 @@ def __eq__(self, other: object) -> bool:
         """Check if two references are equal."""
         if not isinstance(other, Reference):
             return False
-        return type(self) == type(other) and self.vertices == other.vertices
+        return type(self) is type(other) and self.vertices == other.vertices
 
     def __hash__(self) -> int:
         """Check if two references are equal."""

test/utils.py

@@ -22,6 +22,7 @@
         "Wu-Xu": [({}, [2])],
         "MWX": [({}, [1])],
         "EG": [({}, range(1, 5))],
+        "P1-iso-P2": [({}, [1])],
     },
     "triangle": {
         "P": [({"variant": "equispaced"}, range(5))],