From f340caa3fb64f3d348517544999cc67bc3308d8a Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 8 Aug 2023 12:32:29 +0100 Subject: [PATCH] Add P1-iso-P2 on interval (#232) * Add P1 iso P2 on interval * changelog * is not == * piecewisefunction * mypy * 1d piecewise functions * flake * isort --- CHANGELOG_SINCE_LAST_VERSION.md | 1 + README.md | 1 + symfem/elements/p1_iso_p2.py | 54 ++++++++++++++++++++++++++++----- symfem/finite_element.py | 5 ++- symfem/geometry.py | 19 ++++++++++++ symfem/piecewise_functions.py | 7 ++++- symfem/references.py | 2 +- test/utils.py | 1 + 8 files changed, 79 insertions(+), 11 deletions(-) diff --git a/CHANGELOG_SINCE_LAST_VERSION.md b/CHANGELOG_SINCE_LAST_VERSION.md index af1ce586..4935ff1d 100644 --- a/CHANGELOG_SINCE_LAST_VERSION.md +++ b/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 diff --git a/README.md b/README.md index 9629042a..43fd5c90 100644 --- a/README.md +++ b/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) diff --git a/symfem/elements/p1_iso_p2.py b/symfem/elements/p1_iso_p2.py index aea3b277..1ca7194f 100644 --- a/symfem/elements/p1_iso_p2.py +++ b/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 = [] diff --git a/symfem/finite_element.py b/symfem/finite_element.py index dae830ea..833adfce 100644 --- a/symfem/finite_element.py +++ b/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: diff --git a/symfem/geometry.py b/symfem/geometry.py index 49618738..33cb958e 100644 --- a/symfem/geometry.py +++ b/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. diff --git a/symfem/piecewise_functions.py b/symfem/piecewise_functions.py index ac8ddb22..68c2ebaa 100644 --- a/symfem/piecewise_functions.py +++ b/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: diff --git a/symfem/references.py b/symfem/references.py index b498f75d..7e251d9d 100644 --- a/symfem/references.py +++ b/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.""" diff --git a/test/utils.py b/test/utils.py index ffc96845..7f66a33f 100644 --- a/test/utils.py +++ b/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))],