diff --git a/fuse/triples.py b/fuse/triples.py index ce2b412..fe98685 100644 --- a/fuse/triples.py +++ b/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,7 +381,7 @@ 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): @@ -383,7 +389,6 @@ def generate(self, cell, space, id_counter): 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] diff --git a/test/eigen.py b/test/eigen.py index 2fba062..f4ed702 100644 --- a/test/eigen.py +++ b/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) diff --git a/test/element_examples.py b/test/element_examples.py index 4307e30..f3b78a7 100644 --- a/test/element_examples.py +++ b/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") diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index fe6f6f1..2689419 100644 --- a/test/test_convert_to_fiat.py +++ b/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')),