Add orthonormal polynomials (#190)

* Add function to create orthonormal functions

* changelog

Author: mscroggs

CHANGELOG_SINCE_LAST_VERSION.md

@@ -1 +1,2 @@
 - Remove svgwrite dependency
+- Added normalised orthogonal functions

symfem/polynomials.py

@@ -779,3 +779,18 @@ def orthogonal_basis(
         return orthogonal_basis_pyramid(*args)
 
     raise ValueError(f"Unsupported cell type: {cell}")
+
+
+def orthonormal_basis(
+    cell, order: int, derivs: int, variables: AxisVariablesNotSingle = None
+) -> typing.List[typing.List[ScalarFunction]]:
+    """Create a basis of orthonormal polynomials."""
+    from .create import create_reference
+
+    poly = orthogonal_basis(cell, order, derivs, variables)
+    ref = create_reference(cell)
+    norms = [sympy.sqrt((f ** 2).integral(ref, x)) for f in poly[0]]
+    for i, n in enumerate(norms):
+        for j in range(len(poly)):
+            poly[j][i] /= n
+    return poly

test/test_polynomials.py

@@ -7,7 +7,8 @@
 
 from symfem import create_element, create_reference
 from symfem.functions import VectorFunction
-from symfem.polynomials import Hcurl_polynomials, Hdiv_polynomials, orthogonal_basis
+from symfem.polynomials import (Hcurl_polynomials, Hdiv_polynomials, orthogonal_basis,
+                                orthonormal_basis)
 from symfem.symbols import t, x
 
 
@@ -120,3 +121,18 @@ def test_orthogonal_polynomial_derivatives(reference, order):
     for i, j, k in second_d:
         for p, q in zip(polynomials[0], polynomials[i]):
             assert p.diff(x[j]).diff(x[k]) == q
+
+
+@pytest.mark.parametrize("reference", [
+    "interval", "triangle", "quadrilateral",
+    "tetrahedron", "hexahedron", "prism",
+    "pyramid"])
+@pytest.mark.parametrize("order", range(3))
+def test_orthonormal_polynomials(reference, order, speed):
+    if speed == "fast" and order > 2:
+        pytest.skip()
+
+    polynomials = orthonormal_basis(reference, order, 0)[0]
+    ref = create_reference(reference)
+    for p in polynomials:
+        assert (p ** 2).integral(ref, x) == 1