Work on Crouzeix Raviart Elements

fuse/triples.py

@@ -291,7 +291,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set):
                             oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = np.eye(len(ent_dofs_ids))
                         elif g in dof_gen_class.g1.members():
                             sub_mat = g.matrix_form()
-                            oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy()
+                            generators = list(set([ed.sub_id for ed in ent_dofs]))
+                            for gen in generators:
+                                sub_ids = [ed.id for ed in ent_dofs if ed.sub_id == gen]
+                                oriented_mats_by_entity[dim][e_id][val][np.ix_(sub_ids, sub_ids)] = sub_mat.copy()
                         else:
                             pass
                             # # sub component dense case
@@ -328,7 +331,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set):
                     elif len(dof_gen_class.keys()) == 1:
                         if g in dof_gen_class[dim].g1.members() or (pure_perm and len(dof_gen_class[dim].g1.members()) > 1):
                             sub_mat = g.matrix_form()
-                            oriented_mats_overall[val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy()
+                            generators = list(set([ed.sub_id for ed in ent_dofs]))
+                            for gen in generators:
+                                sub_ids = [ed.id for ed in ent_dofs if ed.sub_id == gen]
+                                oriented_mats_overall[val][np.ix_(sub_ids, sub_ids)] = sub_mat.copy()
 
         for val, mat in oriented_mats_overall.items():
             cell_dofs = entity_ids[dim][0]
@@ -375,15 +381,14 @@ def generate(self, cell, space, id_counter):
         if self.ls is None:
             self.ls = []
             for l_g in self.x:
-                i = 0
+                i = id_counter
                 for g in self.g1.members():
                     generated = l_g(g)
                     if not isinstance(generated, list):
                         generated = [generated]
                     for dof in generated:
                         dof.add_context(self, cell, space, g, id_counter, i)
                         id_counter += 1
-                        i += 1
                     self.ls.extend(generated)
             self.dof_numbers = len(self.ls)
             self.dof_ids = [dof.id for dof in self.ls]

test/eigen.py

@@ -3,6 +3,7 @@
 import numpy as np
 from test_2d_examples_docs import construct_cg3
 from test_convert_to_fiat import create_cr, create_cg1
+from element_examples import CR_n, CG_n
 
 
 def create_cr3(cell):
@@ -26,12 +27,14 @@ def create_cr3(cell):
 cr3 = create_cr3(polygon(3))
 cr1 = create_cr(polygon(3))
 cg1 = create_cg1(polygon(3))
+cg5 = CG_n(polygon(3), 5)
+cr5 = CR_n(polygon(3), 5)
 
 for N in [50, 100, 200]:
     mesh = RectangleMesh(N, N, pi, pi)
 
-    for elem, space in zip([cg3, cr3], ["CG", "CR"]):
-        V = FunctionSpace(mesh, elem.to_ufl_elem())
+    for elem, space in zip([cg5, cr5], ["CG", "CR"]):
+        V = FunctionSpace(mesh, elem.to_ufl())
         u = TrialFunction(V)
         v = TestFunction(V)
 

test/element_examples.py

@@ -38,15 +38,19 @@ def convert_to_generation(coords, verts=[(-1, -np.sqrt(3)/3), (0, 2*np.sqrt(3)/3
     divide_1 = ((verts[0][0] + verts[1][0])/2, (verts[0][1] + verts[1][1])/2)
     divide_2 = ((verts[0][0] + verts[2][0])/2, (verts[0][1] + verts[2][1])/2)
     for coord in coords:
-        if check_multiple(coord, verts[0]) or (check_multiple(coord, divide_1) and check_below_line(divide_2, (0, 0), coord) <= 0) or (check_multiple(coord, divide_2) and check_below_line(divide_1, (0, 0), coord) <= 0):
-            coords_C3 += [coord]
-        elif check_below_line(verts[0], (0, 0), coord) == -1 and check_below_line(divide_2, (0, 0), coord) == -1:
+        if check_on_line(verts[0], (0, 0), coord) == -1 and check_on_line(divide_2, (0, 0), coord) == -1:
             coords_S3 += [coord]
+        elif check_multiple(coord, verts[0]) and check_on_line(divide_1, (0, 0), coord) == -1 and check_on_line(divide_2, (0, 0), coord) == -1:
+            coords_C3 += [coord]
+        elif check_multiple(coord, divide_1) and check_on_line(divide_2, (0, 0), coord) == -1:
+            coords_C3 += [coord]
+        elif check_multiple(coord, divide_2) and check_on_line(divide_1, (0, 0), coord) == -1:
+            coords_C3 += [coord]
     assert n == len(coords_S1) + len(coords_S3)*6 + len(coords_C3)*3
     return coords_S1, coords_C3, coords_S3
 
 
-def check_below_line(seg_1, seg_2, coord):
+def check_on_line(seg_1, seg_2, coord):
     if seg_1[0] - seg_2[0] == 0:
         if coord[0] == seg_1[0]:
             return 0
@@ -76,10 +80,12 @@ def check_below_line(seg_1, seg_2, coord):
 
 
 def check_multiple(coord_1, coord_2):
-    return check_below_line(coord_2, (0, 0), coord_1) == 0
+    return check_on_line(coord_2, (0, 0), coord_1) == 0
 
 
 def CR_n(cell, deg):
+    if deg % 2 == 0:
+        raise ValueError("Non-Conforming CR only well defined for odd orders")
     points = np.polynomial.legendre.leggauss(deg)[0]
     Pk = PolynomialSpace(deg)
     sym_points = [DOF(DeltaPairing(), PointKernel((pt,))) for pt in points[:len(points)//2]]
@@ -97,8 +103,43 @@ def CR_n(cell, deg):
     c3 = [DOF(DeltaPairing(), PointKernel(c)) for c in c3]
     s3 = [DOF(DeltaPairing(), PointKernel(c)) for c in s3]
 
-    return ElementTriple(cell, (Pk, CellL2, C0), [DOFGenerator(edge_xs, C3, S1), DOFGenerator(s1, S1, S1), DOFGenerator(c3, C3, S1), DOFGenerator(s3, S3, S1)])
+    generators = [DOFGenerator(edge_xs, C3, S1)]
+    if len(s1) > 0:
+        generators += [DOFGenerator(s1, S1, S1)]
+    if len(c3) > 0:
+        generators += [DOFGenerator(c3, C3, S1)]
+    if len(s3) > 0:
+        generators += [DOFGenerator(s3, S3, S1)]
+    return ElementTriple(cell, (Pk, CellL2, C0), generators)
+
+    # return ElementTriple(cell, (Pk, CellL2, C0), [DOFGenerator(edge_xs, C3, S1), DOFGenerator(s1, S1, S1), DOFGenerator(c3, C3, S1), DOFGenerator(s3, S3, S1)])
+
+
+def CG_n(cell, deg):
+    xs = [DOF(DeltaPairing(), PointKernel(()))]
+    dg0 = ElementTriple(cell.vertices()[0], (P0, CellL2, C0), DOFGenerator(xs, S1, S1))
+
+    v_xs = [immerse(cell, dg0, TrH1)]
+    v_dofs = DOFGenerator(v_xs, C3, S1)
+
+    points = np.linspace(-1, 1, deg + 1)[1:-1]
+    Pk = PolynomialSpace(deg)
+    sym_points = [DOF(DeltaPairing(), PointKernel((pt,))) for pt in points[:len(points)//2]]
+    if 0 in points:
+        edge_dg0 = ElementTriple(cell.edges()[0], (Pk, CellL2, C0), [DOFGenerator(sym_points, S2, S1),
+                                                                     DOFGenerator([DOF(DeltaPairing(), PointKernel((0,)))], S1, S1)])
+    else:
+        edge_dg0 = ElementTriple(cell.edges()[0], (Pk, CellL2, C0), [DOFGenerator(sym_points, S2, S1)])
+    edge_xs = [immerse(cell, edge_dg0, TrH1)]
+
+    interior_coords = triangle_coords(triangle_nums[deg - 3])
+    s1, c3, s3 = convert_to_generation(interior_coords)
+
+    s1 = [DOF(DeltaPairing(), PointKernel(c)) for c in s1]
+    c3 = [DOF(DeltaPairing(), PointKernel(c)) for c in c3]
+    s3 = [DOF(DeltaPairing(), PointKernel(c)) for c in s3]
 
+    return ElementTriple(cell, (Pk, CellL2, C0), [DOFGenerator(edge_xs, C3, S1), DOFGenerator(s1, S1, S1), DOFGenerator(c3, C3, S1), DOFGenerator(s3, S3, S1), v_dofs])
 
 # coords = triangle_coords(21)
 # edge_coords = [(-1, -np.sqrt(3)/3), (0, 2*np.sqrt(3)/3), (1, -np.sqrt(3)/3)]
@@ -115,8 +156,9 @@ def CR_n(cell, deg):
 # ax.scatter([e[0] for e in c3], [e[1] for e in c3], color="green")
 # ax.scatter([e[0] for e in s3], [e[1] for e in s3], color="blue")
 # # ax.scatter([e[0] for e in edge_coords], [e[1] for e in edge_coords], color="blue")
-# # ax.figure.savefig("triangle_gen.png")
-# crn = CR_n(polygon(3), 5)
+# ax.figure.savefig("triangle_gen.png")
+# crn = CG_n(polygon(3), 5)
 # print(len(crn.generate()))
 # for d in crn.generate():
 #     print(d)
+# crn.plot(filename="cg5.png")

test/test_convert_to_fiat.py

@@ -7,7 +7,7 @@
 from test_2d_examples_docs import construct_nd, construct_rt, construct_cg3
 from test_3d_examples_docs import construct_tet_rt
 from test_polynomial_space import flatten
-from element_examples import CR_n
+from element_examples import CR_n, CG_n
 
 
 def create_dg1(cell):
@@ -276,7 +276,13 @@ def test_2d(elem_gen, elem_code, deg):
     assert np.allclose(res, 0)
 
 
-@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [(create_cg1, "CG", 1, 1.8), (create_cg2_tri, "CG", 2, 2.8)])
+@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [(create_cg1, "CG", 1, 1.8),
+                                                              (create_cg2_tri, "CG", 2, 2.8),
+                                                              #   (construct_cg3, "CG", 3, 3.8),
+                                                              #   (lambda cell: CG_n(cell, 3), "CG", 3, 3.8),
+                                                              #   (lambda cell: CG_n(cell, 4), "CG", 4, 4.8)])
+                                                              pytest.param(lambda cell: CG_n(cell, 3), "CG", 3, 3.8, marks=pytest.mark.xfail(reason='Need generic orientations ? or other reason unsure')),
+                                                              pytest.param(lambda cell: CG_n(cell, 4), "CG", 4, 4.8, marks=pytest.mark.xfail(reason='Need generic orientations ? or other reason unsure'))])
 def test_helmholtz(elem_gen, elem_code, deg, conv_rate):
     cell = polygon(3)
     elem = elem_gen(cell)
@@ -459,8 +465,9 @@ def test_quad(params, elem_gen):
                                                     (create_cg1, "CG", 1),
                                                     (create_dg1, "DG", 1),
                                                     (create_cr, "CR", 1),
-                                                    (create_cr3, "CR", 1),
-                                                    (lambda cell: CR_n(cell, 3), "CR", 1),
+                                                    (create_cr3, "CR", 1),  # higher order
+                                                    (lambda cell: CR_n(cell, 5), "CR", 1),
+                                                    (lambda cell: CG_n(cell, 5), "CG", 5),
                                                     (create_cf, "CR", 1),  # Don't think Crouzeix Falk in in Firedrake
                                                     (construct_cg3, "CG", 3),
                                                     pytest.param(construct_nd, "N1curl", 1, marks=pytest.mark.xfail(reason='Dense Matrices needed')),