Merged in reference-geometry (pull request #24)

Minor reference cell improvements

FIAT/reference_element.py

@@ -197,6 +197,13 @@ class Simplex(Cell):
 
     def compute_normal(self, facet_i):
         """Returns the unit normal vector to facet i of codimension 1."""
+        # Interval case
+        if self.get_shape() == LINE:
+            verts = numpy.asarray(self.vertices)
+            v_i, = self.get_topology()[0][facet_i]
+            n = verts[v_i] - verts[[1, 0][v_i]]
+            return n / numpy.linalg.norm(n)
+
         # first, let's compute the span of the simplex
         # This is trivial if we have a d-simplex in R^d.
         # Not so otherwise.
@@ -349,7 +356,7 @@ def make_points(self, dim, entity_id, order):
             return image_pts
 
     def volume(self):
-        """Computes the volumne of the simplex in the appropriate
+        """Computes the volume of the simplex in the appropriate
         dimensional measure."""
         return volume(self.get_vertices())
 
@@ -360,15 +367,16 @@ def volume_of_subcomplex(self, dim, facet_no):
     def compute_scaled_normal(self, facet_i):
         """Returns the unit normal to facet_i of scaled by the
         volume of that facet."""
-        t = self.get_topology()
-        sd = self.get_spatial_dimension()
-        facet_verts_ids = t[sd - 1][facet_i]
-        facet_verts_coords = self.get_vertices_of_subcomplex(facet_verts_ids)
-
-        v = volume(facet_verts_coords)
-
+        dim = self.get_spatial_dimension()
+        v = self.volume_of_subcomplex(dim - 1, facet_i)
         return self.compute_normal(facet_i) * v
 
+    def compute_reference_normal(self, facet_dim, facet_i):
+        """Returns the unit normal in infinity norm to facet_i."""
+        assert facet_dim == self.get_spatial_dimension() - 1
+        n = Simplex.compute_normal(self, facet_i)  # skip UFC overrides
+        return n / numpy.linalg.norm(n, numpy.inf)
+
     def get_facet_transform(self, facet_i):
         """Return a function f such that for a point with facet coordinates
         x_f on facet_i, x_c = f(x_f) is the corresponding cell coordinates.
@@ -714,6 +722,25 @@ def transform(point):
                                 for t, s in zip(sct, slices)]))
         return transform
 
+    def volume(self):
+        """Computes the volume in the appropriate dimensional measure."""
+        return numpy.prod([c.volume() for c in self.cells])
+
+    def compute_reference_normal(self, facet_dim, facet_i):
+        """Returns the unit normal in infinity norm to facet_i of
+        subelement dimension facet_dim."""
+        assert len(facet_dim) == len(self.get_dimension())
+        indicator = numpy.array(self.get_dimension()) - numpy.array(facet_dim)
+        (cell_i,), = numpy.nonzero(indicator)
+
+        n = []
+        for i, c in enumerate(self.cells):
+            if cell_i == i:
+                n.extend(c.compute_reference_normal(facet_dim[i], facet_i))
+            else:
+                n.extend([0] * c.get_spatial_dimension())
+        return numpy.asarray(n)
+
     def contains_point(self, point, epsilon=0):
         """Checks if reference cell contains given point
         (with numerical tolerance)."""
@@ -776,6 +803,19 @@ def get_entity_transform(self, dim, entity_i):
                 2: lambda e: ((1, 1), e)}[dim](entity_i)
         return self.product.get_entity_transform(d, e)
 
+    def volume(self):
+        """Computes the volume in the appropriate dimensional measure."""
+        return self.product.volume()
+
+    def compute_reference_normal(self, facet_dim, facet_i):
+        """Returns the unit normal in infinity norm to facet_i."""
+        assert facet_dim == 1
+        d, i = {0: ((0, 1), 0),
+                1: ((0, 1), 1),
+                2: ((1, 0), 0),
+                3: ((1, 0), 1)}[facet_i]
+        return self.product.compute_reference_normal(d, i)
+
     def contains_point(self, point, epsilon=0):
         """Checks if reference cell contains given point
         (with numerical tolerance)."""

test/unit/test_reference_element.py

@@ -0,0 +1,101 @@
+# Copyright (C) 2016 Miklos Homolya
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+
+from __future__ import absolute_import, print_function, division
+
+import pytest
+import numpy as np
+
+from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron
+from FIAT.reference_element import Point, FiredrakeQuadrilateral, TensorProductCell
+
+point = Point()
+interval = UFCInterval()
+triangle = UFCTriangle()
+quadrilateral = FiredrakeQuadrilateral()
+tetrahedron = UFCTetrahedron()
+interval_x_interval = TensorProductCell(interval, interval)
+triangle_x_interval = TensorProductCell(triangle, interval)
+quadrilateral_x_interval = TensorProductCell(quadrilateral, interval)
+
+
+@pytest.mark.parametrize(('cell', 'volume'),
+                         [pytest.mark.xfail((point, 1)),
+                          (interval, 1),
+                          (triangle, 1/2),
+                          (quadrilateral, 1),
+                          (tetrahedron, 1/6),
+                          (interval_x_interval, 1),
+                          (triangle_x_interval, 1/2),
+                          (quadrilateral_x_interval, 1)])
+def test_volume(cell, volume):
+    assert np.allclose(volume, cell.volume())
+
+
+@pytest.mark.parametrize(('cell', 'normals'),
+                         [(interval, [[-1],
+                                      [1]]),
+                          (triangle, [[1, 1],
+                                      [-1, 0],
+                                      [0, -1]]),
+                          (quadrilateral, [[-1, 0],
+                                           [1, 0],
+                                           [0, -1],
+                                           [0, 1]]),
+                          (tetrahedron, [[1, 1, 1],
+                                         [-1, 0, 0],
+                                         [0, -1, 0],
+                                         [0, 0, -1]])])
+def test_reference_normal(cell, normals):
+    facet_dim = cell.get_spatial_dimension() - 1
+    for facet_number in range(len(cell.get_topology()[facet_dim])):
+        assert np.allclose(normals[facet_number],
+                           cell.compute_reference_normal(facet_dim, facet_number))
+
+
+@pytest.mark.parametrize('cell',
+                         [interval_x_interval,
+                          triangle_x_interval,
+                          quadrilateral_x_interval])
+def test_reference_normal_horiz(cell):
+    dim = cell.get_spatial_dimension()
+    np.allclose((0,) * (dim - 1) + (-1,),
+                cell.compute_reference_normal((dim - 1, 0), 0))  # bottom facet
+    np.allclose((0,) * (dim - 1) + (1,),
+                cell.compute_reference_normal((dim - 1, 0), 1))  # top facet
+
+
+@pytest.mark.parametrize(('cell', 'normals'),
+                         [(interval_x_interval, [[-1, 0],
+                                                 [1, 0]]),
+                          (triangle_x_interval, [[1, 1, 0],
+                                                 [-1, 0, 0],
+                                                 [0, -1, 0]]),
+                          (quadrilateral_x_interval, [[-1, 0, 0],
+                                                      [1, 0, 0],
+                                                      [0, -1, 0],
+                                                      [0, 1, 0]])])
+def test_reference_normal_vert(cell, normals):
+    dim = cell.get_spatial_dimension()
+    vert_dim = (dim - 2, 1)
+    for facet_number in range(len(cell.get_topology()[vert_dim])):
+        assert np.allclose(normals[facet_number],
+                           cell.compute_reference_normal(vert_dim, facet_number))
+
+if __name__ == '__main__':
+    import os
+    pytest.main(os.path.abspath(__file__))