add SymmetricSimplex

FIAT/reference_element.py

@@ -687,6 +687,36 @@ def distance_to_point_l1(self, point):
         return abs(l1_dist)
 
 
+class DefaultSimplex(Simplex):
+
+    def get_facet_element(self):
+        dimension = self.get_spatial_dimension()
+        return self.construct_subelement(dimension - 1)
+
+    def construct_subelement(self, dimension):
+        """Constructs the reference element of a cell subentity
+        specified by subelement dimension.
+
+        :arg dimension: subentity dimension (integer)
+        """
+        return default_simplex(dimension)
+
+
+class SymmetricSimplex(Simplex):
+
+    def get_facet_element(self):
+        dimension = self.get_spatial_dimension()
+        return self.construct_subelement(dimension - 1)
+
+    def construct_subelement(self, dimension):
+        """Constructs the reference element of a cell subentity
+        specified by subelement dimension.
+
+        :arg dimension: subentity dimension (integer)
+        """
+        return symmetric_simplex(dimension)
+
+
 class Point(Simplex):
     """This is the reference point."""
 
@@ -705,7 +735,7 @@ def construct_subelement(self, dimension):
         return self
 
 
-class DefaultLine(Simplex):
+class DefaultLine(DefaultSimplex):
     """This is the reference line with vertices (-1.0,) and (1.0,)."""
 
     def __init__(self):
@@ -715,9 +745,6 @@ def __init__(self):
                     1: edges}
         super(DefaultLine, self).__init__(LINE, verts, topology)
 
-    def get_facet_element(self):
-        raise NotImplementedError()
-
 
 class UFCInterval(UFCSimplex):
     """This is the reference interval with vertices (0.0,) and (1.0,)."""
@@ -730,7 +757,7 @@ def __init__(self):
         super(UFCInterval, self).__init__(LINE, verts, topology)
 
 
-class DefaultTriangle(Simplex):
+class DefaultTriangle(DefaultSimplex):
     """This is the reference triangle with vertices (-1.0,-1.0),
     (1.0,-1.0), and (-1.0,1.0)."""
 
@@ -744,9 +771,6 @@ def __init__(self):
                     1: edges, 2: faces}
         super(DefaultTriangle, self).__init__(TRIANGLE, verts, topology)
 
-    def get_facet_element(self):
-        return DefaultLine()
-
 
 class UFCTriangle(UFCSimplex):
     """This is the reference triangle with vertices (0.0,0.0),
@@ -786,7 +810,7 @@ def get_facet_element(self):
         return UFCInterval()
 
 
-class DefaultTetrahedron(Simplex):
+class DefaultTetrahedron(DefaultSimplex):
     """This is the reference tetrahedron with vertices (-1,-1,-1),
     (1,-1,-1),(-1,1,-1), and (-1,-1,1)."""
 
@@ -811,9 +835,6 @@ def __init__(self):
         topology = {0: vs, 1: edges, 2: faces, 3: tets}
         super(DefaultTetrahedron, self).__init__(TETRAHEDRON, verts, topology)
 
-    def get_facet_element(self):
-        return DefaultTriangle()
-
 
 class IntrepidTetrahedron(Simplex):
     """This is the reference tetrahedron with vertices (0,0,0),
@@ -1308,6 +1329,18 @@ def ufc_simplex(spatial_dim):
         raise RuntimeError("Can't create UFC simplex of dimension %s." % str(spatial_dim))
 
 
+def symmetric_simplex(spatial_dim):
+    A = numpy.array([[2, 1, 1],
+                     [0, numpy.sqrt(3), numpy.sqrt(3)/3],
+                     [0, 0, numpy.sqrt(6)*(2/3)]])
+    A = A[:spatial_dim, :][:, :spatial_dim]
+    b = A.sum(axis=1) * (-1 / (1 + spatial_dim))
+    Ref1 = ufc_simplex(spatial_dim)
+    v = numpy.dot(Ref1.get_vertices(), A.T) + b[None, :]
+    vertices = tuple(map(tuple, v))
+    return SymmetricSimplex(Ref1.get_shape(), vertices, Ref1.get_topology())
+
+
 def ufc_cell(cell):
     """Handle incoming calls from FFC."""