From efb644cab410fb4f7c58337cce2ec5ca08bf5096 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 7 Jan 2025 17:12:32 +0000 Subject: [PATCH 1/2] sorta working general cg and cr, tbc on convergence rates --- redefining_fe/triples.py | 13 +++++--- test/eigen.py | 7 +++-- test/element_examples.py | 58 +++++++++++++++++++++++++++++++----- test/test_convert_to_fiat.py | 9 +++--- 4 files changed, 69 insertions(+), 18 deletions(-) diff --git a/redefining_fe/triples.py b/redefining_fe/triples.py index 7e1378e..37f60d1 100644 --- a/redefining_fe/triples.py +++ b/redefining_fe/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 696ecd3..c2be0d7 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 1bca580..8c7122b 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,7 +80,7 @@ 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): @@ -97,8 +101,45 @@ 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): + if deg % 2 == 0: + raise ValueError("Non-Conforming CR only well defined for odd orders") + 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 ff757b6..e36c7a6 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,7 @@ 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), pytest.param(lambda cell: CG_n(cell, 3), "CG", 3, 3.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 +459,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')), From f6f875687e53c16a323c715a3660d54523ca6930 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 13 Jan 2025 12:26:03 +0000 Subject: [PATCH 2/2] Error in where error was being raised --- test/element_examples.py | 4 ++-- test/test_convert_to_fiat.py | 8 +++++++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/test/element_examples.py b/test/element_examples.py index d7d9baf..f3b78a7 100644 --- a/test/element_examples.py +++ b/test/element_examples.py @@ -84,6 +84,8 @@ def check_multiple(coord_1, coord_2): 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]] @@ -114,8 +116,6 @@ def CR_n(cell, deg): def CG_n(cell, deg): - if deg % 2 == 0: - raise ValueError("Non-Conforming CR only well defined for odd orders") xs = [DOF(DeltaPairing(), PointKernel(()))] dg0 = ElementTriple(cell.vertices()[0], (P0, CellL2, C0), DOFGenerator(xs, S1, S1)) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index f700470..2689419 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -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.param(lambda cell: CG_n(cell, 3), "CG", 3, 3.8, marks=pytest.mark.xfail(reason='Need generic orientations ? or other reason unsure')),]) +@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)