change reference cell

having a different reference cell will strengthen the documentation

Author: mscroggs

README.md

@@ -32,7 +32,7 @@ values = e.tabulate(points)
 
 ## Conventions
 ### Reference cell
-The reference cell used by simplefem is the triangle with vertices at (0, 0), (1, 0),
+The reference cell used by simplefem is the triangle with vertices at (-1, 0), (1, 0),
 and (0, 1).
 
 ### Point ordering
@@ -41,13 +41,13 @@ are ordered lexicographically: for example, the points that define a degree 3 el
 are arranged like this:
 
 ```
-9
-|\
-7 8
-|  \
-4 5 6
-|    \
-0-1-2-3
+      9
+     / \
+    7   8
+   /     \
+  4   5   6
+ /         \
+0---1---2---3
 ```
 
 

simplefem/lagrange.py

@@ -12,7 +12,11 @@ def __init__(self, degree):
         """Initialise."""
         self.degree = degree
         self.evaluation_points = np.array(
-            [[j / degree, i / degree] for i in range(degree + 1) for j in range(degree + 1 - i)]
+            [
+                [(2 * j + i) / degree - 1, i / degree]
+                for i in range(degree + 1)
+                for j in range(degree + 1 - i)
+            ]
         )
 
         self.coeffs = np.linalg.inv(tabulate(self.evaluation_points, degree))

simplefem/polynomials.py

@@ -44,7 +44,7 @@ def tabulate(points: npt.NDArray[np.float64], degree: int) -> npt.NDArray[np.flo
     table = np.zeros(((degree + 1) * (degree + 2) // 2, points.shape[0]))
     table[_index(0, 0)] = 1
 
-    x = points[:, 0]
+    x = (points[:, 0] - points[:, 1] + 1) / 2
     y = points[:, 1]
 
     for p in range(1, degree + 1):

test/test_lagrange.py

@@ -13,8 +13,8 @@ def test_lagrange1_element():
     x = points[:, 0]
     y = points[:, 1]
 
-    assert np.allclose(values[0], 1 - x - y)
-    assert np.allclose(values[1], x)
+    assert np.allclose(values[0], (1 - x - y) / 2)
+    assert np.allclose(values[1], (1 + x - y) / 2)
     assert np.allclose(values[2], y)
 
 
@@ -27,9 +27,16 @@ def test_lagrange2_element():
     x = points[:, 0]
     y = points[:, 1]
 
-    assert np.allclose(values[0], (1 - x - y) * (1 - 2 * x - 2 * y))
-    assert np.allclose(values[1], 4 * x * (1 - x - y))
-    assert np.allclose(values[2], x * (2 * x - 1))
-    assert np.allclose(values[3], 4 * y * (1 - x - y))
-    assert np.allclose(values[4], 4 * x * y)
+    assert np.allclose(values[0], (y + x - 1) * (y + x) / 2)
+    assert np.allclose(values[1], (1 + x - y) * (1 - x - y))
+    assert np.allclose(values[2], (y - x - 1) * (y - x) / 2)
+    assert np.allclose(values[3], 2 * y * (1 - x - y))
+    assert np.allclose(values[4], 2 * y * (1 + x - y))
     assert np.allclose(values[5], y * (2 * y - 1))
+
+
+def test_vertices():
+    e = lagrange_element(1)
+    assert np.allclose(e.evaluation_points[0], [-1, 0])
+    assert np.allclose(e.evaluation_points[1], [1, 0])
+    assert np.allclose(e.evaluation_points[2], [0, 1])

test/test_orthogonal.py

@@ -11,7 +11,7 @@ def test_orthogonal(degree):
         quadraturerules.Domain.Triangle,
         max(1, 2 * degree),
     )
-    pts = pts[:, 1:]
+    pts = np.array([[p[1] - p[0], p[2]] for p in pts])
     values = tabulate(pts, degree)
 
     for i, v_i in enumerate(values):
@@ -26,7 +26,7 @@ def test_orthonormal(degree):
         quadraturerules.Domain.Triangle,
         max(1, 2 * degree),
     )
-    pts = pts[:, 1:]
+    pts = np.array([[p[1] - p[0], p[2]] for p in pts])
     values = tabulate(pts, degree)
 
     for v in values:
@@ -35,7 +35,7 @@ def test_orthonormal(degree):
 
 @pytest.mark.parametrize("degree", range(6))
 def test_not_nan_at_vertices(degree):
-    pts = np.array([[0, 0], [1, 0], [0, 1]])
+    pts = np.array([[-1, 0], [1, 0], [0, 1]])
 
     values = tabulate(pts, degree)