From 80304b7d7ebe8d8b5962ea1dee56eb93b6d60ecb Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 1 Apr 2025 14:09:27 +0100 Subject: [PATCH 01/32] add equality to cells --- fuse/__init__.py | 3 ++- fuse/cells.py | 22 +++++++++++++++++++++- test/test_cells.py | 10 ++++++++++ test/test_convert_to_fiat.py | 6 ++++-- test/test_tensor_prod.py | 12 ++++++++++++ 5 files changed, 49 insertions(+), 4 deletions(-) diff --git a/fuse/__init__.py b/fuse/__init__.py index e7ed0618..9b68648d 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,4 +1,5 @@ -from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex + +from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex, TensorProductPoint from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, diff_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, VectorKernel, PolynomialKernel, ComponentKernel from fuse.triples import ElementTriple, DOFGenerator, immerse diff --git a/fuse/cells.py b/fuse/cells.py index ab06eb98..340b1035 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -151,6 +151,13 @@ def compute_scaled_verts(d, n): raise ValueError("Dimension {} not supported".format(d)) +def line(): + """ + Constructs the default 1D interval + """ + return Point(1, [Point(0), Point(0)], vertex_num=2) + + def polygon(n): """ Constructs the 2D default cell with n sides/vertices @@ -544,6 +551,9 @@ def ordered_vertices(self, get_class=False): return self.oriented.permute(verts) return verts + def ordered_vertex_coords(self): + return [self.get_node(o, return_coords=True) for o in self.ordered_vertices()] + def d_entities_ids(self, d): return self.d_entities(d, get_class=False) @@ -645,7 +655,6 @@ def basis_vectors(self, return_coords=True, entity=None, order=False): self_levels = [generation for generation in nx.topological_generations(self.G)] vertices = entity.ordered_vertices() if self.dimension == 0: - # return [[] raise ValueError("Dimension 0 entities cannot have Basis Vectors") if self.oriented: # ordered_vertices() handles the orientation so we want to drop the orientation node @@ -874,6 +883,16 @@ def dict_id(self): def _from_dict(o_dict): return Point(o_dict["dim"], o_dict["edges"], oriented=o_dict["oriented"], cell_id=o_dict["id"]) + def __eq__(self, other): + if self.dimension != other.dimension: + return False + if set(self.ordered_vertex_coords()) != set(other.ordered_vertex_coords()): + return False + return self.get_topology() == other.get_topology() + + def __hash__(self): + return hash(self.id) + class Edge(): """ @@ -971,6 +990,7 @@ def to_fiat(self, name=None): return CellComplexToFiatTensorProduct(self, name) def flatten(self): + assert self.A == self.B return TensorProductPoint(self.A, self.B, True) diff --git a/test/test_cells.py b/test/test_cells.py index ced8e36d..54b9aa5a 100644 --- a/test/test_cells.py +++ b/test/test_cells.py @@ -179,6 +179,16 @@ def test_comparison(): # print(tensor_product1 >= tensor_product1) +def test_self_equality(C): + assert C == C + + +@pytest.mark.parametrize(["A", "B", "res"], [(firedrake_triangle(), polygon(3), False), + (line(), line(), True),]) +def test_equality(A, B, res): + assert (A == B) == res + + @pytest.mark.parametrize(["cell"], [(ufc_triangle(),), (polygon(3),)]) def test_connectivity(cell): cell = cell.to_fiat() diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index bee86500..af195fb0 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -542,8 +542,10 @@ def test_quad(elem_gen): assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) -def test_non_tensor_quad(): - create_cg1_quad() +# def test_non_tensor_quad(): +# elem = create_cg1_quad() +# ufl_elem = elem.to_ufl() +# assert (run_test(0, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) def project(U, mesh, func): diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 7c527203..67232963 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -117,3 +117,15 @@ def test_quad_mesh_helmholtz(): conv = np.log2(res[:-1] / res[1:]) print("convergence order:", conv) assert (np.array(conv) > 1.8).all() + + +@pytest.mark.parametrize(["A", "B", "res"], [(Point(0), line(), False), + (line(), line(), True), + (polygon(3), line(), False),]) +def test_flattening(A, B, res): + tensor_cell = TensorProductPoint(A, B) + if not res: + with pytest.raises(AssertionError): + tensor_cell.flatten() + else: + tensor_cell.flatten() From d881c4a89f9bd696726d5f04cea5dd0411470d9f Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 1 Apr 2025 16:16:58 +0100 Subject: [PATCH 02/32] refactor --- fuse/cells.py | 29 ++++++++++++++++++++++------- test/test_convert_to_fiat.py | 8 ++++---- 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 340b1035..10d850d5 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -948,11 +948,11 @@ def _from_dict(o_dict): class TensorProductPoint(): - def __init__(self, A, B, flat=False): + def __init__(self, A, B): self.A = A self.B = B self.dimension = self.A.dimension + self.B.dimension - self.flat = flat + self.flat = False def ordered_vertices(self): return self.A.ordered_vertices() + self.B.ordered_vertices() @@ -980,18 +980,33 @@ def vertices(self, get_class=True, return_coords=False): return verts def to_ufl(self, name=None): - if self.flat: - return CellComplexToUFL(self, "quadrilateral") return TensorProductCell(self.A.to_ufl(), self.B.to_ufl()) def to_fiat(self, name=None): - if self.flat: - return CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name)) return CellComplexToFiatTensorProduct(self, name) def flatten(self): assert self.A == self.B - return TensorProductPoint(self.A, self.B, True) + return FlattenedPoint(self.A, self.B) + + +class FlattenedPoint(TensorProductPoint): + + def __init__(self, A, B): + self.A = A + self.B = B + self.dimension = self.A.dimension + self.B.dimension + self.flat = True + + def to_ufl(self, name=None): + return CellComplexToUFL(self, "quadrilateral") + + def to_fiat(self, name=None): + # TODO this should check if it actually is a hypercube + return CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name)) + + def flatten(self): + return self class CellComplexToFiatSimplex(Simplex): diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index af195fb0..598bb9dd 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -542,10 +542,10 @@ def test_quad(elem_gen): assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) -# def test_non_tensor_quad(): -# elem = create_cg1_quad() -# ufl_elem = elem.to_ufl() -# assert (run_test(0, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) +def test_non_tensor_quad(): + elem = create_cg1_quad() + ufl_elem = elem.to_ufl() + assert (run_test(0, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) def project(U, mesh, func): From e49c7034544b06c7921b7b235cfc9b75b6896075 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 1 Apr 2025 17:10:00 +0100 Subject: [PATCH 03/32] first steps to making a fuse element from a flattened tensor prod --- fuse/cells.py | 24 ++++++++++++++++++++++-- test/test_tensor_prod.py | 3 ++- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 10d850d5..bc43781a 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -964,8 +964,8 @@ def get_sub_entities(self): self.A.get_sub_entities() self.B.get_sub_entities() - def dimension(self): - return tuple(self.A.dimension, self.B.dimension) + def dim(self): + return (self.A.dimension, self.B.dimension) def d_entities(self, d, get_class=True): return self.A.d_entities(d, get_class) + self.B.d_entities(d, get_class) @@ -1005,6 +1005,26 @@ def to_fiat(self, name=None): # TODO this should check if it actually is a hypercube return CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name)) + def construct_fuse_rep(self): + sub_cells = [self.A, self.B] + dims = self.dim() + points = {i: [] for i in range(max(dims))} + + for d in range(max(dims)): + points = [] + attachments = [] + for cell in sub_cells: + if d <= cell.dimension: + points.append(cell.d_entities(d)) + attachments.append(cell.edges()) + for i in itertools.product(*points): + print(i) + print(attachments) + new_attachments = [] + # for a in attachments: + # old_attach = a.attachment + # new_attachments.append() + def flatten(self): return self diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 67232963..d8b2bba8 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -128,4 +128,5 @@ def test_flattening(A, B, res): with pytest.raises(AssertionError): tensor_cell.flatten() else: - tensor_cell.flatten() + cell = tensor_cell.flatten() + cell.construct_fuse_rep() \ No newline at end of file From 2fb6b88abd5ffff5bee0813b79733f2c6ff53037 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 3 Apr 2025 14:14:08 +0100 Subject: [PATCH 04/32] first level of generated hasse diagram --- fuse/cells.py | 94 ++++++++++++++++++++++++++++------- test/test_cells.py | 120 +-------------------------------------------- 2 files changed, 79 insertions(+), 135 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index bc43781a..c7f5a536 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -509,6 +509,15 @@ def get_starter_ids(self): min_ids = [min(dimension) for dimension in structure] return min_ids + def local_id(self, node): + structure = [sorted(generation) for generation in nx.topological_generations(self.G)] + structure.reverse() + min_id = self.get_starter_ids() + for d in range(len(structure)): + if node.id in structure[d]: + return node.id - min_id[d] + raise ValueError("Node not found in cell") + def graph_dim(self): if self.oriented: dim = self.dimension + 1 @@ -882,16 +891,14 @@ def dict_id(self): def _from_dict(o_dict): return Point(o_dict["dim"], o_dict["edges"], oriented=o_dict["oriented"], cell_id=o_dict["id"]) - - def __eq__(self, other): + + def equivalent(self, other): if self.dimension != other.dimension: return False if set(self.ordered_vertex_coords()) != set(other.ordered_vertex_coords()): return False return self.get_topology() == other.get_topology() - def __hash__(self): - return hash(self.id) class Edge(): @@ -986,7 +993,7 @@ def to_fiat(self, name=None): return CellComplexToFiatTensorProduct(self, name) def flatten(self): - assert self.A == self.B + assert self.A.equivalent(self.B) return FlattenedPoint(self.A, self.B) @@ -1008,22 +1015,75 @@ def to_fiat(self, name=None): def construct_fuse_rep(self): sub_cells = [self.A, self.B] dims = self.dim() - points = {i: [] for i in range(max(dims))} + points = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells} + attachments = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells} - for d in range(max(dims)): - points = [] - attachments = [] + for d in range(max(dims) + 1): for cell in sub_cells: + + # if d <= cell.dimension: + # points[d].extend([cell.get_node(e, return_coords=True) for e in cell.d_entities(d, get_class=False)]) if d <= cell.dimension: - points.append(cell.d_entities(d)) - attachments.append(cell.edges()) - for i in itertools.product(*points): - print(i) + sub_ent = cell.d_entities(d, get_class=True) + + # points[d].extend([cell.get_node(e, return_coords=True) for e in cell.d_entities(d, get_class=False)]) + points[cell][d].extend(sub_ent) + for s in sub_ent: + attachments[cell][d].extend(s.connections) + + new_attachments = [] + for a in attachments[cell][d]: + print(cell) + print(a, " res: ", a.attachment) print(attachments) - new_attachments = [] - # for a in attachments: - # old_attach = a.attachment - # new_attachments.append() + print(points) + + + # old_attach = a.attachment + # new_attachments.append() + prod_points = list(itertools.product(*[points[cell][0] for cell in sub_cells])) + print(len(prod_points)) + point_cls = [Point(0) for i in range(len(prod_points))] + edges = [] + + # generate edges of tensor product result + for a in prod_points: + for b in prod_points: + # of all combinations of point, take those where at least one changes and at least one is the same + if any(a[i] == b[i] for i in range(len(a))) and any(a[i] != b[i] for i in range(len(sub_cells))): + # ensure if they change, that edge exists in the existing topology + if all([a[i]== b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i].topology[1].values()) for i in range(len(sub_cells))]): + edges.append((a, b)) + + # hasse level 1 + edge_cls1 = {e: [] for e in edges} + print(prod_points) + for i in range(len(sub_cells)): + for (a, b) in edges: + a_idx = prod_points.index(a) + b_idx = prod_points.index(b) + if a[i] != b[i]: + a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] + b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0] + edge_cls1[(a,b)].append(Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), + Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)])) + print(edge_cls1) + # hasse level 2 + # for i in range(len(sub_cells)): + # for (a, b) in edges: + # a_idx = prod_points.index(a) + # b_idx = prod_points.index(b) + # print(a, a_idx) + # print(b, b_idx) + # # breakpoint() + # if a[i] != b[i]: + # a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] + # b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0] + # edge_cls1.append(Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), + # Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)])) + # print(attachments) + # for i in itertools.product(attachments, repeat=2): + # print(i) def flatten(self): return self diff --git a/test/test_cells.py b/test/test_cells.py index 54b9aa5a..cac57755 100644 --- a/test/test_cells.py +++ b/test/test_cells.py @@ -185,121 +185,5 @@ def test_self_equality(C): @pytest.mark.parametrize(["A", "B", "res"], [(firedrake_triangle(), polygon(3), False), (line(), line(), True),]) -def test_equality(A, B, res): - assert (A == B) == res - - -@pytest.mark.parametrize(["cell"], [(ufc_triangle(),), (polygon(3),)]) -def test_connectivity(cell): - cell = cell.to_fiat() - for dim0 in range(cell.get_spatial_dimension()+1): - connectivity = cell.get_connectivity()[(dim0, 0)] - topology = cell.get_topology()[dim0] - assert len(connectivity) == len(topology) - - assert all(connectivity[i] == t for i, t in topology.items()) - - -def test_tensor_connectivity(): - from test_2d_examples_docs import construct_cg1 - A = construct_cg1() - B = construct_cg1() - cell = tensor_product(A, B).cell - cell = cell.to_fiat() - for dim0 in [(0, 0), (1, 0), (0, 1), (1, 1)]: - connectivity = cell.get_connectivity()[(dim0, (0, 0))] - topology = cell.get_topology()[dim0] - assert len(connectivity) == len(topology) - - assert all(connectivity[i] == t for i, t in topology.items()) - - -@pytest.mark.parametrize(["cell"], [(ufc_triangle(),), (polygon(3),), (make_tetrahedron(), ), (make_tetrahedron(), )]) -def test_new_connectivity(cell): - cell = cell.to_fiat() - for dim0 in range(cell.get_dimension() + 1): - connectivity = cell.get_connectivity()[(dim0, 0)] - topology = cell.get_topology()[dim0] - assert len(connectivity) == len(topology) - for i, t in topology.items(): - print(connectivity[i]) - print(t) - assert all(connectivity[i] == t for i, t in topology.items()) - - -def test_compare_tris(): - fuse_tet = polygon(3) - ufc_tet = ufc_triangle() - fiat_tet = ufc_simplex(2) - - print(fiat_tet.get_topology()) - print(fuse_tet.get_topology()) - print(ufc_tet.get_topology()) - fiat_connectivity = fiat_tet.get_connectivity() - fuse_connectivity = fuse_tet.to_fiat().get_connectivity() - ufc_connectivity = ufc_tet.to_fiat().get_connectivity() - _dim = fiat_tet.get_dimension() - print("fiat") - print(make_entity_cone_lists(fiat_tet)) - for dim0 in range(_dim): - connectivity = fiat_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse") - print(make_entity_cone_lists(fuse_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = fuse_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse ufc") - print(make_entity_cone_lists(ufc_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = ufc_connectivity[(dim0+1, dim0)] - print(connectivity) - - -def test_compare_tets(): - tet = make_tetrahedron() - # perm = tet.group.get_member([1, 2, 0, 3]) - fuse_tet = tet - ufc_tet = ufc_tetrahedron() - fiat_tet = ufc_simplex(3) - # breakpoint() - print(fiat_tet.get_topology()) - print(fuse_tet.get_topology()) - print(ufc_tet.get_topology()) - fiat_connectivity = fiat_tet.get_connectivity() - fuse_connectivity = fuse_tet.to_fiat().get_connectivity() - ufc_connectivity = ufc_tet.to_fiat().get_connectivity() - _dim = fiat_tet.get_dimension() - print("fiat") - print(make_entity_cone_lists(fiat_tet)) - for dim0 in range(_dim): - connectivity = fiat_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse") - print(make_entity_cone_lists(fuse_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = fuse_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse ufc") - print(make_entity_cone_lists(ufc_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = ufc_connectivity[(dim0+1, dim0)] - print(connectivity) - - -def make_entity_cone_lists(fiat_cell): - _dim = fiat_cell.get_dimension() - _connectivity = fiat_cell.connectivity - _list = [] - _offset_list = [0 for _ in _connectivity[(0, 0)]] # vertices have no cones - _offset = 0 - _n = 0 # num. of entities up to dimension = _d - for _d in range(_dim): - _n1 = len(_offset_list) - for _conn in _connectivity[(_d + 1, _d)]: - _list += [_c + _n for _c in _conn] # These are indices into cell_closure[some_cell] - _offset_list.append(_offset) - _offset += len(_conn) - _n = _n1 - _offset_list.append(_offset) - return _list, _offset_list +def test_equivalence(A, B, res): + assert A.equivalent(B) == res From e6341b9e7a0996dd493233471bda2991d8ab92ba Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 7 Apr 2025 12:25:50 +0100 Subject: [PATCH 05/32] fix recursion issue --- fuse/cells.py | 56 +++++++++++++++++++----------------- test/test_convert_to_fiat.py | 3 +- 2 files changed, 32 insertions(+), 27 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index c7f5a536..732b41ce 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -380,6 +380,7 @@ def compute_cell_group(self): """ verts = self.ordered_vertices() v_coords = [self.get_node(v, return_coords=True) for v in verts] + n = len(verts) max_group = SymmetricGroup(n) edges = [edge.ordered_vertices() for edge in self.edges()] @@ -923,7 +924,11 @@ def __call__(self, *x): if hasattr(self.attachment, '__iter__'): res = [] for attach_comp in self.attachment: - res.append(sympy_to_numpy(attach_comp, syms, x)) + if len(attach_comp.atoms(sp.Symbol)) <= len(x): + res.append(sympy_to_numpy(attach_comp, syms, x)) + else: + res.append(attach_comp.subs({syms[i]: x[i] for i in range(len(x))})) + return tuple(res) return sympy_to_numpy(self.attachment, syms, x) return x @@ -997,13 +1002,15 @@ def flatten(self): return FlattenedPoint(self.A, self.B) -class FlattenedPoint(TensorProductPoint): +class FlattenedPoint(Point, TensorProductPoint): def __init__(self, A, B): self.A = A self.B = B self.dimension = self.A.dimension + self.B.dimension self.flat = True + fuse_edges = self.construct_fuse_rep() + super().__init__(self.dimension, fuse_edges) def to_ufl(self, name=None): return CellComplexToUFL(self, "quadrilateral") @@ -1014,7 +1021,8 @@ def to_fiat(self, name=None): def construct_fuse_rep(self): sub_cells = [self.A, self.B] - dims = self.dim() + dims = (self.A.dimension, self.B.dimension) + points = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells} attachments = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells} @@ -1038,7 +1046,6 @@ def construct_fuse_rep(self): print(attachments) print(points) - # old_attach = a.attachment # new_attachments.append() prod_points = list(itertools.product(*[points[cell][0] for cell in sub_cells])) @@ -1056,7 +1063,7 @@ def construct_fuse_rep(self): edges.append((a, b)) # hasse level 1 - edge_cls1 = {e: [] for e in edges} + edge_cls1 = {e: None for e in edges} print(prod_points) for i in range(len(sub_cells)): for (a, b) in edges: @@ -1065,25 +1072,23 @@ def construct_fuse_rep(self): if a[i] != b[i]: a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0] - edge_cls1[(a,b)].append(Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), - Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)])) - print(edge_cls1) + edge_cls1[(a,b)] = Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), + Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)]) + edge_cls2 = [] # hasse level 2 - # for i in range(len(sub_cells)): - # for (a, b) in edges: - # a_idx = prod_points.index(a) - # b_idx = prod_points.index(b) - # print(a, a_idx) - # print(b, b_idx) - # # breakpoint() - # if a[i] != b[i]: - # a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] - # b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0] - # edge_cls1.append(Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), - # Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)])) - # print(attachments) - # for i in itertools.product(attachments, repeat=2): - # print(i) + for i in range(len(sub_cells)): + for (a, b) in edges: + if a[i] == b[i]: + x = sp.Symbol("x") + a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] + if i == 0: + attach = (x,) + a_edge.attachment + else: + attach = a_edge.attachment + (x,) + print(edge_cls1[(a, b)].ordered_vertices()) + edge_cls2.append(Edge(edge_cls1[(a, b)], attach, a_edge.o)) + print(edge_cls2) + return edge_cls2 def flatten(self): return self @@ -1276,9 +1281,8 @@ def constructCellComplex(name): return polygon(3).to_ufl(name) # return ufc_triangle().to_ufl(name) elif name == "quadrilateral": - interval = Point(1, [Point(0), Point(0)], vertex_num=2) - return TensorProductPoint(interval, interval).flatten().to_ufl(name) - # return ufc_quad().to_ufl(name) + return TensorProductPoint(line(), line()).flatten().to_ufl(name) + # return firedrake_quad().to_ufl(name) # return polygon(4).to_ufl(name) elif name == "tetrahedron": # return ufc_tetrahedron().to_ufl(name) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 598bb9dd..3caf5ee2 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -117,6 +117,7 @@ def create_cg1_quad(): xs = [immerse(cell, vert_dg, TrH1)] Pk = PolynomialSpace(deg, deg + 1) + breakpoint() cg = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, get_cyc_group(len(cell.vertices())), S1)) return cg @@ -545,7 +546,7 @@ def test_quad(elem_gen): def test_non_tensor_quad(): elem = create_cg1_quad() ufl_elem = elem.to_ufl() - assert (run_test(0, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) + assert (run_test(1, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) def project(U, mesh, func): From 208d64a77fb98cc47db3969d1be571f17bee9c39 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 7 Apr 2025 12:27:17 +0100 Subject: [PATCH 06/32] lint --- fuse/cells.py | 15 ++++++--------- test/test_tensor_prod.py | 2 +- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 732b41ce..e0bed57b 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -380,7 +380,7 @@ def compute_cell_group(self): """ verts = self.ordered_vertices() v_coords = [self.get_node(v, return_coords=True) for v in verts] - + n = len(verts) max_group = SymmetricGroup(n) edges = [edge.ordered_vertices() for edge in self.edges()] @@ -892,7 +892,7 @@ def dict_id(self): def _from_dict(o_dict): return Point(o_dict["dim"], o_dict["edges"], oriented=o_dict["oriented"], cell_id=o_dict["id"]) - + def equivalent(self, other): if self.dimension != other.dimension: return False @@ -901,7 +901,6 @@ def equivalent(self, other): return self.get_topology() == other.get_topology() - class Edge(): """ Representation of the connections in a cell complex. @@ -928,7 +927,6 @@ def __call__(self, *x): res.append(sympy_to_numpy(attach_comp, syms, x)) else: res.append(attach_comp.subs({syms[i]: x[i] for i in range(len(x))})) - return tuple(res) return sympy_to_numpy(self.attachment, syms, x) return x @@ -1039,7 +1037,6 @@ def construct_fuse_rep(self): for s in sub_ent: attachments[cell][d].extend(s.connections) - new_attachments = [] for a in attachments[cell][d]: print(cell) print(a, " res: ", a.attachment) @@ -1059,9 +1056,9 @@ def construct_fuse_rep(self): # of all combinations of point, take those where at least one changes and at least one is the same if any(a[i] == b[i] for i in range(len(a))) and any(a[i] != b[i] for i in range(len(sub_cells))): # ensure if they change, that edge exists in the existing topology - if all([a[i]== b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i].topology[1].values()) for i in range(len(sub_cells))]): + if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i].topology[1].values()) for i in range(len(sub_cells))]): edges.append((a, b)) - + # hasse level 1 edge_cls1 = {e: None for e in edges} print(prod_points) @@ -1072,8 +1069,8 @@ def construct_fuse_rep(self): if a[i] != b[i]: a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0] - edge_cls1[(a,b)] = Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), - Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)]) + edge_cls1[(a, b)] = Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), + Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)]) edge_cls2 = [] # hasse level 2 for i in range(len(sub_cells)): diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index d8b2bba8..8fa82f27 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -129,4 +129,4 @@ def test_flattening(A, B, res): tensor_cell.flatten() else: cell = tensor_cell.flatten() - cell.construct_fuse_rep() \ No newline at end of file + cell.construct_fuse_rep() From 3f06f030f3ed9509d2c826e87a0ae539fb46fd22 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 7 Apr 2025 17:59:55 +0100 Subject: [PATCH 07/32] refactor work on poly spaces and trying to figure out how to make quad spaces --- fuse/spaces/polynomial_spaces.py | 57 ++++++++++++++++++++++++-------- fuse/utils.py | 21 +++++++----- test/test_convert_to_fiat.py | 20 +++++++++-- test/test_polynomial_space.py | 12 +++++++ 4 files changed, 86 insertions(+), 24 deletions(-) diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 422cf053..0b74783b 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -3,7 +3,7 @@ from FIAT.reference_element import cell_to_simplex from FIAT import expansions, polynomial_set, reference_element from itertools import chain -from fuse.utils import tabulate_sympy, max_deg_sp_mat +from fuse.utils import tabulate_sympy, max_deg_sp_expr import sympy as sp import numpy as np from functools import total_ordering @@ -56,9 +56,27 @@ def to_ON_polynomial_set(self, ref_el, k=None): shape = (sd,) else: shape = tuple() + base_ON = ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal") + + # if self.contains: + # expansion_set = expansions.ExpansionSet(ref_el) + # num_exp_functions = expansion_set.get_num_members(self.maxdegree) + # dimPmin = expansions.polynomial_dimension(ref_el, self.mindegree) + # dimPmax = expansions.polynomial_dimension(ref_el, self.maxdegree) + # indices = list(chain(*(range(i * dimPmin, i * dimPmax) for i in range(sd)))) + # from math import comb + # print(list(chain(*[range(0, 0), range(2, 3)]))) + # base_ON.expansion_set._tabulate_on_cell(3, np.array([[0,0]])) + # if self.set_shape: + # indices = list(chain(*(range(i * dimPmin, i * dimPmax) for i in range(sd)))) + # else: + # indices = list(range(dimPmin, dimPmax)) + # restricted_ON = base_ON.take(indices) + # return restricted_ON + + # breakpoint() if self.mindegree > 0: - base_ON = ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal") dimPmin = expansions.polynomial_dimension(ref_el, self.mindegree) dimPmax = expansions.polynomial_dimension(ref_el, self.maxdegree) if self.set_shape: @@ -67,7 +85,8 @@ def to_ON_polynomial_set(self, ref_el, k=None): indices = list(range(dimPmin, dimPmax)) restricted_ON = base_ON.take(indices) return restricted_ON - return ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal") + + return base_ON def __repr__(self): res = "" @@ -163,37 +182,49 @@ def to_ON_polynomial_set(self, ref_el): if not isinstance(ref_el, reference_element.Cell): ref_el = ref_el.to_fiat() k = max([s.maxdegree for s in self.spaces]) - space_poly_sets = [s.to_ON_polynomial_set(ref_el) for s in self.spaces] sd = ref_el.get_spatial_dimension() ref_el = cell_to_simplex(ref_el) - if all([w == 1 for w in self.weights]): - weighted_sets = space_poly_sets - # otherwise have to work on this through tabulation Q = create_quadrature(ref_el, 2 * (k + 1)) Qpts, Qwts = Q.get_points(), Q.get_weights() weighted_sets = [] - for (space, w) in zip(space_poly_sets, self.weights): + for (s, w) in zip(self.spaces, self.weights): + space = s.to_ON_polynomial_set(ref_el) + if s.set_shape: + shape = (sd,) + else: + shape = tuple() if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)): weighted_sets.append(space) else: - w_deg = max_deg_sp_mat(w) - Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, scale="orthonormal") - vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal") + if isinstance(w, sp.Expr): + w = sp.Matrix([[w]]) + vec = False + else: + vec = True + w_deg = max_deg_sp_expr(w) + Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, shape, scale="orthonormal") + # vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal") space_at_Qpts = space.tabulate(Qpts)[(0,) * sd] Pkpw_at_Qpts = Pkpw.tabulate(Qpts)[(0,) * sd] tabulated_expr = tabulate_sympy(w, Qpts).T - scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] + if s.set_shape or vec: + scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] + else: + # breakpoint() + scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] + scaled_at_Qpts = scaled_at_Qpts.squeeze() PkHw_coeffs = np.dot(np.multiply(scaled_at_Qpts, Qwts), Pkpw_at_Qpts.T) + # breakpoint() weighted_sets.append(polynomial_set.PolynomialSet(ref_el, space.degree + w_deg, space.degree + w_deg, - vec_Pkpw.get_expansion_set(), + Pkpw.get_expansion_set(), PkHw_coeffs)) combined_sets = weighted_sets[0] for i in range(1, len(weighted_sets)): diff --git a/fuse/utils.py b/fuse/utils.py index fcfa3d9e..3a91c6ba 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -47,7 +47,7 @@ def tabulate_sympy(expr, pts): # expr: sp matrix expression in x,y,z for components of R^d # pts: n values in R^d # returns: evaluation of expr at pts - res = np.array(pts) + res = np.zeros((pts.shape[0],) + (expr.shape[-1],)) i = 0 syms = ["x", "y", "z"] for pt in pts: @@ -57,16 +57,21 @@ def tabulate_sympy(expr, pts): subbed = np.array(subbed).astype(np.float64) res[i] = subbed[0] i += 1 - final = res.squeeze() - return final + # final = res.squeeze() + return res -def max_deg_sp_mat(sp_mat): +def max_deg_sp_expr(sp_expr): degs = [] - for comp in sp_mat: - # only compute degree if component is a polynomial - if sp.sympify(comp).as_poly(): - degs += [sp.sympify(comp).as_poly().degree()] + if isinstance(sp_expr, sp.Matrix): + for comp in sp_expr: + # only compute degree if component is a polynomial + if sp.sympify(comp).as_poly(): + degs += [sp.sympify(comp).as_poly().degree()] + else: + if sp.sympify(sp_expr).as_poly(): + degs += [sp.sympify(sp_expr).as_poly().degree()] + return max(degs) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 3caf5ee2..cdf0dd5d 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -115,9 +115,9 @@ def create_cg1_quad(): vert_dg = create_dg0(cell.vertices()[0]) xs = [immerse(cell, vert_dg, TrH1)] + x = sp.Symbol("y") - Pk = PolynomialSpace(deg, deg + 1) - breakpoint() + Pk = PolynomialSpace(deg) + P2.restrict(1, 1)*x cg = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, get_cyc_group(len(cell.vertices())), S1)) return cg @@ -535,7 +535,7 @@ def test_poisson_analytic(params, elem_gen): @pytest.mark.parametrize(['elem_gen'], - [(create_cg1_quad_tensor,), pytest.param(create_cg1_quad, marks=pytest.mark.xfail(reason='Need to allow generation on tensor product quads'))]) + [(create_cg1_quad_tensor,), pytest.param(create_cg1_quad, marks=pytest.mark.xfail(reason='How to make quad poly set correctly'))]) def test_quad(elem_gen): elem = elem_gen() r = 0 @@ -563,6 +563,20 @@ def project(U, mesh, func): return res +def project(U, mesh, func): + f = assemble(interpolate(func, U)) + + out = Function(U) + u = TrialFunction(U) + v = TestFunction(U) + a = inner(u, v)*dx + L = inner(f, v)*dx + solve(a == L, out) + + res = sqrt(assemble(dot(out - func, out - func) * dx)) + return res + + @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), (create_cg1, "CG", 1), (create_dg1, "DG", 1), diff --git a/test/test_polynomial_space.py b/test/test_polynomial_space.py index e7dad66d..d22ebd54 100644 --- a/test/test_polynomial_space.py +++ b/test/test_polynomial_space.py @@ -41,6 +41,7 @@ def test_restriction(): res_on_set = restricted.to_ON_polynomial_set(cell) P3_on_set = P3.to_ON_polynomial_set(cell) + assert res_on_set.get_num_members() < P3_on_set.get_num_members() not_restricted = P3.restrict(0, 3) @@ -48,6 +49,17 @@ def test_restriction(): assert not_restricted.mindegree == 0 +@pytest.mark.xfail(reason="Quad space WIP") +def test_square_space(): + cell = polygon(3) + q2 = PolynomialSpace(3, 1) + + q2_on_set = q2.to_ON_polynomial_set(cell) + P3_on_set = P3.to_ON_polynomial_set(cell) + + assert q2_on_set.get_num_members() < P3_on_set.get_num_members() + + @pytest.mark.parametrize("deg", [1, 2, 3, 4]) def test_complete_space(deg): cell = polygon(3) From 01ac541c93c68278c4a18abb0e8bd114db04f94a Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 10 Apr 2025 12:19:33 +0100 Subject: [PATCH 08/32] fix contains in polynomial sets --- fuse/spaces/polynomial_spaces.py | 34 +++++++++++--------------------- test/test_convert_to_fiat.py | 34 ++------------------------------ test/test_polynomial_space.py | 1 - 3 files changed, 14 insertions(+), 55 deletions(-) diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 0b74783b..3d4618fd 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -1,4 +1,5 @@ from FIAT.polynomial_set import ONPolynomialSet +from FIAT.expansions import morton_index2, morton_index3 from FIAT.quadrature_schemes import create_quadrature from FIAT.reference_element import cell_to_simplex from FIAT import expansions, polynomial_set, reference_element @@ -8,6 +9,8 @@ import numpy as np from functools import total_ordering +morton_index = {2: morton_index2, 3: morton_index3} + @total_ordering class PolynomialSpace(object): @@ -47,7 +50,6 @@ def degree(self): return self.maxdegree def to_ON_polynomial_set(self, ref_el, k=None): - # how does super/sub degrees work here if not isinstance(ref_el, reference_element.Cell): ref_el = ref_el.to_fiat() sd = ref_el.get_spatial_dimension() @@ -57,24 +59,7 @@ def to_ON_polynomial_set(self, ref_el, k=None): else: shape = tuple() base_ON = ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal") - - # if self.contains: - # expansion_set = expansions.ExpansionSet(ref_el) - # num_exp_functions = expansion_set.get_num_members(self.maxdegree) - # dimPmin = expansions.polynomial_dimension(ref_el, self.mindegree) - # dimPmax = expansions.polynomial_dimension(ref_el, self.maxdegree) - # indices = list(chain(*(range(i * dimPmin, i * dimPmax) for i in range(sd)))) - # from math import comb - # print(list(chain(*[range(0, 0), range(2, 3)]))) - # base_ON.expansion_set._tabulate_on_cell(3, np.array([[0,0]])) - # if self.set_shape: - # indices = list(chain(*(range(i * dimPmin, i * dimPmax) for i in range(sd)))) - # else: - # indices = list(range(dimPmin, dimPmax)) - # restricted_ON = base_ON.take(indices) - # return restricted_ON - - # breakpoint() + indices = None if self.mindegree > 0: dimPmin = expansions.polynomial_dimension(ref_el, self.mindegree) @@ -83,10 +68,15 @@ def to_ON_polynomial_set(self, ref_el, k=None): indices = list(chain(*(range(i * dimPmin, i * dimPmax) for i in range(sd)))) else: indices = list(range(dimPmin, dimPmax)) - restricted_ON = base_ON.take(indices) - return restricted_ON - return base_ON + if self.contains != self.maxdegree and self.contains != -1: + indices = [morton_index[sd](p, q) for p in range(self.contains + 1) for q in range(self.contains + 1)] + + if indices is None: + return base_ON + + restricted_ON = base_ON.take(indices) + return restricted_ON def __repr__(self): res = "" diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index cdf0dd5d..d0a5f641 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -1,6 +1,5 @@ import pytest import numpy as np -import sympy as sp from fuse import * from firedrake import * from sympy.combinatorics import Permutation @@ -115,9 +114,7 @@ def create_cg1_quad(): vert_dg = create_dg0(cell.vertices()[0]) xs = [immerse(cell, vert_dg, TrH1)] - x = sp.Symbol("y") - - Pk = PolynomialSpace(deg) + P2.restrict(1, 1)*x + Pk = PolynomialSpace(deg + 1, deg) cg = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, get_cyc_group(len(cell.vertices())), S1)) return cg @@ -543,40 +540,13 @@ def test_quad(elem_gen): assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) +@pytest.mark.xfail(reason="Issue with quad cell") def test_non_tensor_quad(): elem = create_cg1_quad() ufl_elem = elem.to_ufl() assert (run_test(1, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) -def project(U, mesh, func): - f = assemble(interpolate(func, U)) - - out = Function(U) - u = TrialFunction(U) - v = TestFunction(U) - a = inner(u, v)*dx - L = inner(f, v)*dx - solve(a == L, out) - - res = sqrt(assemble(dot(out - func, out - func) * dx)) - return res - - -def project(U, mesh, func): - f = assemble(interpolate(func, U)) - - out = Function(U) - u = TrialFunction(U) - v = TestFunction(U) - a = inner(u, v)*dx - L = inner(f, v)*dx - solve(a == L, out) - - res = sqrt(assemble(dot(out - func, out - func) * dx)) - return res - - @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), (create_cg1, "CG", 1), (create_dg1, "DG", 1), diff --git a/test/test_polynomial_space.py b/test/test_polynomial_space.py index d22ebd54..baf81206 100644 --- a/test/test_polynomial_space.py +++ b/test/test_polynomial_space.py @@ -49,7 +49,6 @@ def test_restriction(): assert not_restricted.mindegree == 0 -@pytest.mark.xfail(reason="Quad space WIP") def test_square_space(): cell = polygon(3) q2 = PolynomialSpace(3, 1) From 78e0861ddcf8b4588816e5609b23c70f04ed7ae7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 14 Apr 2025 17:20:58 +0100 Subject: [PATCH 09/32] working on cell --- fuse/cells.py | 24 +++++------------------- test/test_convert_to_fiat.py | 11 +++++------ 2 files changed, 10 insertions(+), 25 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index e0bed57b..5b4c436b 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -1015,7 +1015,8 @@ def to_ufl(self, name=None): def to_fiat(self, name=None): # TODO this should check if it actually is a hypercube - return CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name)) + fiat = CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name)) + return fiat def construct_fuse_rep(self): sub_cells = [self.A, self.B] @@ -1026,27 +1027,16 @@ def construct_fuse_rep(self): for d in range(max(dims) + 1): for cell in sub_cells: - - # if d <= cell.dimension: - # points[d].extend([cell.get_node(e, return_coords=True) for e in cell.d_entities(d, get_class=False)]) if d <= cell.dimension: sub_ent = cell.d_entities(d, get_class=True) - - # points[d].extend([cell.get_node(e, return_coords=True) for e in cell.d_entities(d, get_class=False)]) points[cell][d].extend(sub_ent) for s in sub_ent: attachments[cell][d].extend(s.connections) - for a in attachments[cell][d]: - print(cell) - print(a, " res: ", a.attachment) - print(attachments) - print(points) - - # old_attach = a.attachment - # new_attachments.append() prod_points = list(itertools.product(*[points[cell][0] for cell in sub_cells])) - print(len(prod_points)) + # temp = prod_points[1] + # prod_points[1] = prod_points[2] + # prod_points[2] = temp point_cls = [Point(0) for i in range(len(prod_points))] edges = [] @@ -1058,10 +1048,8 @@ def construct_fuse_rep(self): # ensure if they change, that edge exists in the existing topology if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i].topology[1].values()) for i in range(len(sub_cells))]): edges.append((a, b)) - # hasse level 1 edge_cls1 = {e: None for e in edges} - print(prod_points) for i in range(len(sub_cells)): for (a, b) in edges: a_idx = prod_points.index(a) @@ -1082,9 +1070,7 @@ def construct_fuse_rep(self): attach = (x,) + a_edge.attachment else: attach = a_edge.attachment + (x,) - print(edge_cls1[(a, b)].ordered_vertices()) edge_cls2.append(Edge(edge_cls1[(a, b)], attach, a_edge.o)) - print(edge_cls2) return edge_cls2 def flatten(self): diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index d0a5f641..07316d47 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -109,10 +109,9 @@ def create_cg1(cell): def create_cg1_quad(): deg = 1 - cell = polygon(4) - # cell = constructCellComplex("quadrilateral").cell_complex - - vert_dg = create_dg0(cell.vertices()[0]) + cell = TensorProductPoint(line(), line()).flatten() + print(cell, type(cell)) + vert_dg = create_dg1(cell.vertices()[0]) xs = [immerse(cell, vert_dg, TrH1)] Pk = PolynomialSpace(deg + 1, deg) cg = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, get_cyc_group(len(cell.vertices())), S1)) @@ -532,7 +531,7 @@ def test_poisson_analytic(params, elem_gen): @pytest.mark.parametrize(['elem_gen'], - [(create_cg1_quad_tensor,), pytest.param(create_cg1_quad, marks=pytest.mark.xfail(reason='How to make quad poly set correctly'))]) + [(create_cg1_quad_tensor,), pytest.param(create_cg1_quad, marks=pytest.mark.xfail(reason='Issue with cell/mesh'))]) def test_quad(elem_gen): elem = elem_gen() r = 0 @@ -540,7 +539,7 @@ def test_quad(elem_gen): assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) -@pytest.mark.xfail(reason="Issue with quad cell") +# @pytest.mark.xfail(reason="Issue with quad cell") def test_non_tensor_quad(): elem = create_cg1_quad() ufl_elem = elem.to_ufl() From e9560f9c4d227516526568e26255569b17443291 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 23 Apr 2025 12:34:51 +0100 Subject: [PATCH 10/32] Ensure all numeric reps are 0..n over the whole group --- fuse/groups.py | 18 +++++++++++++++++- test/test_convert_to_fiat.py | 13 ++++++++----- test/test_groups.py | 7 +++++++ 3 files changed, 32 insertions(+), 6 deletions(-) diff --git a/fuse/groups.py b/fuse/groups.py index e63a9d32..1ec5275c 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -65,9 +65,15 @@ def compute_perm(self, base_val=None): return val, val_list def numeric_rep(self): + """ Uses a standard formula to number permutations in the group. + For the case where this doesn't automatically number from 0..n (ie the group is not the full symmetry group), + a mapping is constructed on group creation""" identity = self.group.identity.perm.array_form m_array = self.perm.array_form - return orientation_value(identity, m_array) + val = orientation_value(identity, m_array) + if self.group.group_rep_numbering is not None: + return self.group.group_rep_numbering[val] + return val def __eq__(self, x): assert isinstance(x, GroupMemberRep) @@ -144,6 +150,11 @@ def __init__(self, perm_list, cell=None): counter += 1 # self._members = sorted(self._members, key=lambda g: g.numeric_rep()) + self.group_rep_numbering = None + numeric_reps = [m.numeric_rep() for m in self.members()] + if sorted(numeric_reps) != list(range(len(numeric_reps))): + self.group_rep_numbering = {a: b for a, b in zip(sorted(numeric_reps), list(range(len(numeric_reps))))} + def add_cell(self, cell): return PermutationSetRepresentation(self.perm_list, cell=cell) @@ -243,6 +254,11 @@ def __init__(self, base_group, cell=None): self.identity = p_rep counter += 1 + self.group_rep_numbering = None + numeric_reps = [m.numeric_rep() for m in self.members()] + if sorted(numeric_reps) != list(range(len(numeric_reps))): + self.group_rep_numbering = {a: b for a, b in zip(sorted(numeric_reps), list(range(len(numeric_reps))))} + # this order produces simpler generator lists # self.generators.reverse() diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 07316d47..efca5daa 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -500,6 +500,7 @@ def helmholtz_solve(V, mesh): def poisson_solve(r, elem, parameters={}, quadrilateral=False): # Create mesh and define function space m = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) + x = SpatialCoordinate(m) V = FunctionSpace(m, elem) @@ -539,11 +540,13 @@ def test_quad(elem_gen): assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) -# @pytest.mark.xfail(reason="Issue with quad cell") -def test_non_tensor_quad(): - elem = create_cg1_quad() - ufl_elem = elem.to_ufl() - assert (run_test(1, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) +# # @pytest.mark.xfail(reason="Issue with quad cell") +# def test_non_tensor_quad(): +# elem = create_cg1_quad() +# ufl_elem = elem.to_ufl() +# print(elem.to_fiat().entity_permutations()) +# # elem.cell.hasse_diagram(filename="cg1quad.png") +# assert (run_test(1, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), diff --git a/test/test_groups.py b/test/test_groups.py index 1b26787e..0c207039 100644 --- a/test/test_groups.py +++ b/test/test_groups.py @@ -117,3 +117,10 @@ def test_perm_mat_conversion(): mat_form = g.matrix_form() array_form = perm_matrix_to_perm_array(mat_form) assert np.allclose(g.perm.array_form, array_form) + + +def test_numeric_reps(): + cell = polygon(4) + rot4 = get_cyc_group(4).add_cell(cell) + + assert sorted([m.numeric_rep() for m in rot4.members()]) == list(range(len(rot4.members()))) From 37b6124c629399abb9eca29b4291340527e3ff76 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 24 Apr 2025 09:56:52 +0100 Subject: [PATCH 11/32] comparision to existing --- test/test_convert_to_fiat.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index efca5daa..7eb7a2e3 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -521,6 +521,29 @@ def poisson_solve(r, elem, parameters={}, quadrilateral=False): return sqrt(assemble(inner(u - f, u - f) * dx)) +def run_test_original(r, elem_code, deg, parameters={}, quadrilateral=False): + # Create mesh and define function space + m = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) + + x = SpatialCoordinate(m) + V = FunctionSpace(m, elem_code, deg) + # Define variational problem + u = Function(V) + v = TestFunction(V) + a = inner(grad(u), grad(v)) * dx + + bcs = [DirichletBC(V, Constant(0), 3), + DirichletBC(V, Constant(42), 4)] + + # Compute solution + solve(a == 0, u, solver_parameters=parameters, bcs=bcs) + + f = Function(V) + f.interpolate(42*x[1]) + + return sqrt(assemble(inner(u - f, u - f) * dx)) + + @pytest.mark.parametrize(['params', 'elem_gen'], [(p, d) for p in [{}, {'snes_type': 'ksponly', 'ksp_type': 'preonly', 'pc_type': 'lu'}] @@ -546,7 +569,7 @@ def test_quad(elem_gen): # ufl_elem = elem.to_ufl() # print(elem.to_fiat().entity_permutations()) # # elem.cell.hasse_diagram(filename="cg1quad.png") -# assert (run_test(1, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) +# assert (run_test_original(1, "CG", 1, parameters={}, quadrilateral=True) < 1.e-9) @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), From d2d2d4bb18960795ea1b1b1bfc77b743f53489d4 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 10 Mar 2026 18:06:49 +0000 Subject: [PATCH 12/32] tidy up --- fuse/__init__.py | 2 +- fuse/cells.py | 6 ++++-- fuse/spaces/polynomial_spaces.py | 2 +- fuse/utils.py | 2 +- test/test_convert_to_fiat.py | 28 +++++++++++++++++++++------- 5 files changed, 28 insertions(+), 12 deletions(-) diff --git a/fuse/__init__.py b/fuse/__init__.py index 9b68648d..0db1c3d7 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,5 +1,5 @@ -from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex, TensorProductPoint +from fuse.cells import Point, Edge, polygon, line, make_tetrahedron, constructCellComplex, TensorProductPoint from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, diff_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, VectorKernel, PolynomialKernel, ComponentKernel from fuse.triples import ElementTriple, DOFGenerator, immerse diff --git a/fuse/cells.py b/fuse/cells.py index 5b4c436b..b7b25833 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -926,7 +926,9 @@ def __call__(self, *x): if len(attach_comp.atoms(sp.Symbol)) <= len(x): res.append(sympy_to_numpy(attach_comp, syms, x)) else: - res.append(attach_comp.subs({syms[i]: x[i] for i in range(len(x))})) + res_val = attach_comp.subs({syms[i]: x[i] for i in range(len(x))}) + res.append(res_val) + return tuple(res) return sympy_to_numpy(self.attachment, syms, x) return x @@ -1046,7 +1048,7 @@ def construct_fuse_rep(self): # of all combinations of point, take those where at least one changes and at least one is the same if any(a[i] == b[i] for i in range(len(a))) and any(a[i] != b[i] for i in range(len(sub_cells))): # ensure if they change, that edge exists in the existing topology - if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i].topology[1].values()) for i in range(len(sub_cells))]): + if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i]._topology[1].values()) for i in range(len(sub_cells))]): edges.append((a, b)) # hasse level 1 edge_cls1 = {e: None for e in edges} diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 3d4618fd..3502a6ea 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -157,7 +157,7 @@ def __init__(self, weights, spaces): self.weights = weights self.spaces = spaces - weight_degrees = [0 if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)) else max_deg_sp_mat(w) for w in self.weights] + weight_degrees = [0 if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)) else max_deg_sp_expr(w) for w in self.weights] maxdegree = max([space.maxdegree + w_deg for space, w_deg in zip(spaces, weight_degrees)]) mindegree = min([space.mindegree + w_deg for space, w_deg in zip(spaces, weight_degrees)]) diff --git a/fuse/utils.py b/fuse/utils.py index 3a91c6ba..6592e117 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -29,7 +29,7 @@ def sympy_to_numpy(array, symbols, values): """ substituted = array.subs({symbols[i]: values[i] for i in range(len(values))}) - if len(array.atoms(sp.Symbol)) == len(values) and all(not isinstance(v, sp.Expr) for v in values): + if len(array.atoms(sp.Symbol)) <= len(values) and all(not isinstance(v, sp.Expr) for v in values): nparray = np.array(substituted).astype(np.float64) if len(nparray.shape) > 1: diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 7eb7a2e3..2cd8885d 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -563,13 +563,27 @@ def test_quad(elem_gen): assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) -# # @pytest.mark.xfail(reason="Issue with quad cell") -# def test_non_tensor_quad(): -# elem = create_cg1_quad() -# ufl_elem = elem.to_ufl() -# print(elem.to_fiat().entity_permutations()) -# # elem.cell.hasse_diagram(filename="cg1quad.png") -# assert (run_test_original(1, "CG", 1, parameters={}, quadrilateral=True) < 1.e-9) +@pytest.mark.xfail(reason="Issue with quad cell") +def test_non_tensor_quad(): + elem = create_cg1_quad() + ufl_elem = elem.to_ufl() + print(elem.to_fiat().entity_permutations()) + # elem.cell.hasse_diagram(filename="cg1quad.png") + assert (run_test_original(1, "CG", 1, parameters={}, quadrilateral=True) < 1.e-9) + + +def project(U, mesh, func): + f = assemble(interpolate(func, U)) + + out = Function(U) + u = TrialFunction(U) + v = TestFunction(U) + a = inner(u, v)*dx + L = inner(f, v)*dx + solve(a == L, out) + + res = sqrt(assemble(dot(out - func, out - func) * dx)) + return res @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), From 06a50d0f489267085c81c941480292a121f58f7a Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 10 Mar 2026 18:08:11 +0000 Subject: [PATCH 13/32] lint --- test/test_cells.py | 4 ++-- test/test_convert_to_fiat.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_cells.py b/test/test_cells.py index cac57755..0fdbe0e1 100644 --- a/test/test_cells.py +++ b/test/test_cells.py @@ -1,9 +1,9 @@ from fuse import * from firedrake import * -from fuse.cells import ufc_triangle, ufc_tetrahedron +from fuse.cells import ufc_triangle import pytest import numpy as np -from FIAT.reference_element import default_simplex, ufc_simplex +from FIAT.reference_element import default_simplex from test_convert_to_fiat import helmholtz_solve diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 2cd8885d..f0004b32 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -566,7 +566,7 @@ def test_quad(elem_gen): @pytest.mark.xfail(reason="Issue with quad cell") def test_non_tensor_quad(): elem = create_cg1_quad() - ufl_elem = elem.to_ufl() + # ufl_elem = elem.to_ufl() print(elem.to_fiat().entity_permutations()) # elem.cell.hasse_diagram(filename="cg1quad.png") assert (run_test_original(1, "CG", 1, parameters={}, quadrilateral=True) < 1.e-9) From 52fa438d12c2106e21ab024fb8a3661a08204e19 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 10 Mar 2026 18:58:38 +0000 Subject: [PATCH 14/32] small fixes --- test/test_cells.py | 2 +- test/test_convert_to_fiat.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_cells.py b/test/test_cells.py index 0fdbe0e1..c4da9859 100644 --- a/test/test_cells.py +++ b/test/test_cells.py @@ -183,7 +183,7 @@ def test_self_equality(C): assert C == C -@pytest.mark.parametrize(["A", "B", "res"], [(firedrake_triangle(), polygon(3), False), +@pytest.mark.parametrize(["A", "B", "res"], [(ufc_triangle(), polygon(3), False), (line(), line(), True),]) def test_equivalence(A, B, res): assert A.equivalent(B) == res diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index f0004b32..6213fd1d 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -1,5 +1,6 @@ import pytest import numpy as np +import sympy as sp from fuse import * from firedrake import * from sympy.combinatorics import Permutation From 0e197b45e7f987044467a0c32e25c71e71aa81f5 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 10:36:48 +0000 Subject: [PATCH 15/32] modifications to polynomial kernel to take into account value shape --- fuse/dof.py | 37 +++++++++++++++++++++-------------- test/test_2d_examples_docs.py | 32 +++++++++++++++++++++++++++--- test/test_convert_to_fiat.py | 4 +++- 3 files changed, 54 insertions(+), 19 deletions(-) diff --git a/fuse/dof.py b/fuse/dof.py index 7c07f1c9..c576ee48 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -152,7 +152,7 @@ def permute(self, g): def __call__(self, *args): return self.pt - def evaluate(self, Qpts, Qwts, basis_change, immersed, dim): + def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): return np.array([self.pt for _ in Qpts]).astype(np.float64), np.ones_like(Qwts), [[tuple()] for pt in Qpts] def _to_dict(self): @@ -189,12 +189,13 @@ def permute(self, g): def __call__(self, *args): return self.pt - def evaluate(self, Qpts, Qwts, basis_change, immersed, dim): + def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): + comps = [(i,) for v in value_shape for i in range(v)] if isinstance(self.pt, int): - return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] + return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), comps if not immersed: - return Qpts, np.array([wt*np.matmul(self.pt, basis_change)for wt in Qwts]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] - return Qpts, np.array([wt*immersed(np.matmul(self.pt, basis_change))for wt in Qwts]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] + return Qpts, np.array([wt*np.matmul(self.pt, basis_change)for wt in Qwts]).astype(np.float64), comps + return Qpts, np.array([wt*immersed(np.matmul(self.pt, basis_change))for wt in Qwts]).astype(np.float64), comps def _to_dict(self): o_dict = {"pt": self.pt} @@ -209,7 +210,11 @@ def _from_dict(obj_dict): class PolynomialKernel(BaseKernel): - def __init__(self, fn, g=None, symbols=[], shape=0): + def __init__(self, fn, g=None, symbols=[]): + if hasattr(fn, "__iter__"): + shape = len(fn) + else: + shape = 0 if len(symbols) != 0 and (shape != 0 and any(not sp.sympify(fn[i]).as_poly() for i in range(shape))) and not sp.sympify(fn).as_poly(): raise ValueError("Function argument or its components must be able to be interpreted as a sympy polynomial") if shape != 0: @@ -235,19 +240,21 @@ def degree(self, interpolant_degree): def permute(self, g): # new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) new_fn = self.fn - return PolynomialKernel(new_fn, g=g, symbols=self.syms, shape=self.shape) + return PolynomialKernel(new_fn, g=g, symbols=self.syms) def __call__(self, *args): if self.shape == 0: res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) else: - res = [] - for i in range(self.shape): - res += [sympy_to_numpy(self.fn[i], self.syms, args[:len(self.syms)])] + res = [sympy_to_numpy(self.fn[i], self.syms, args[:len(self.syms)]) for i in range(self.shape)] return res - def evaluate(self, Qpts, Qwts, basis_change, immersed, dim): - return Qpts, np.array([wt*self(*(np.matmul(pt, basis_change))) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] + def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): + if len(value_shape) == 0: + comps = [[tuple()] for pt in Qpts] + else: + comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] + return Qpts, np.array([wt*self(*(np.matmul(pt, basis_change))) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), comps def _to_dict(self): o_dict = {"fn": self.fn} @@ -343,9 +350,9 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): # TODO deriv dict needs implementing (currently {}) - return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree), {}, str(self)) + return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree, value_shape), {}, str(self)) - def to_quadrature(self, arg_degree): + def to_quadrature(self, arg_degree, value_shape): Qpts, Qwts = self.cell_defined_on.quadrature(self.kernel.degree(arg_degree)) Qwts = Qwts.reshape(Qwts.shape + (1,)) dim = self.cell_defined_on.get_spatial_dimension() @@ -363,7 +370,7 @@ def immersed(pt): return np.matmul(basis_coeffs, immersed_basis) else: immersed = self.immersed - pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change, immersed, self.cell.dimension) + pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change, immersed, self.cell.dimension, value_shape) if self.immersed: # need to compute jacobian from attachment. diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index a7117416..06da1eea 100644 --- a/test/test_2d_examples_docs.py +++ b/test/test_2d_examples_docs.py @@ -30,6 +30,32 @@ def construct_dg1(): return dg1 +def construct_dg0_integral(cell=None): + edge = Point(1, [Point(0), Point(0)], vertex_num=2) + xs = [DOF(L2Pairing(), PolynomialKernel(1))] + dg0 = ElementTriple(edge, (P0, CellL2, C0), DOFGenerator(xs, S1, S1)) + return dg0 + + +def construct_dg1_integral(cell=None): + edge = Point(1, [Point(0), Point(0)], vertex_num=2) + x = sp.Symbol("x") + xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] + dg1 = ElementTriple(edge, (P1, CellL2, C0), DOFGenerator(xs, S2, S1)) + return dg1 + + +def construct_dg2_integral(cell=None): + edge = Point(1, [Point(0), Point(0)], vertex_num=2) + x = sp.Symbol("x") + xs = [DOF(L2Pairing(), PolynomialKernel((x/2)*(x + 1), symbols=(x,)))] + centre = [DOF(L2Pairing(), PolynomialKernel((1 - x**2), symbols=(x,)))] + + dofs = [DOFGenerator(xs, S2, S1), DOFGenerator(centre, S1, S1)] + dg2 = ElementTriple(edge, (PolynomialSpace(2), CellL2, C0), dofs) + return dg2 + + def plot_dg1(): dg1 = construct_dg1() dg1.plot() @@ -243,9 +269,9 @@ def construct_bdm2(tri=None): phi_0 = [-1/6 - (np.sqrt(3)/6)*y, (-np.sqrt(3)/6) + (np.sqrt(3)/6)*x] phi_1 = [-1/6 - (np.sqrt(3)/6)*y, (np.sqrt(3)/6) + (np.sqrt(3)/6)*x] phi_2 = [1/3 - (np.sqrt(3)/6)*y, (np.sqrt(3)/6)*x] - xs = [DOF(L2Pairing(), PolynomialKernel(phi_0, symbols=(x, y), shape=2)), - DOF(L2Pairing(), PolynomialKernel(phi_1, symbols=(x, y), shape=2)), - DOF(L2Pairing(), PolynomialKernel(phi_2, symbols=(x, y), shape=2))] + xs = [DOF(L2Pairing(), PolynomialKernel(phi_0, symbols=(x, y))), + DOF(L2Pairing(), PolynomialKernel(phi_1, symbols=(x, y))), + DOF(L2Pairing(), PolynomialKernel(phi_2, symbols=(x, y)))] interior = DOFGenerator(xs, S1, S1) nd = PolynomialSpace(deg, set_shape=True) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 6213fd1d..9bc7d3dd 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -5,7 +5,7 @@ from firedrake import * from sympy.combinatorics import Permutation from FIAT.quadrature_schemes import create_quadrature -from test_2d_examples_docs import construct_cg1, construct_nd, construct_rt, construct_cg3 +from test_2d_examples_docs import construct_cg1, construct_nd, construct_rt, construct_cg3, construct_dg1_integral, construct_dg2_integral from test_3d_examples_docs import construct_tet_rt, construct_tet_rt2, construct_tet_ned, construct_tet_ned_2nd_kind, construct_tet_bdm, construct_tet_ned2, construct_tet_cg4 from test_polynomial_space import flatten from element_examples import CR_n @@ -358,6 +358,8 @@ def test_entity_perms(elem_gen, cell): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg1, "CG", 1), (create_dg1, "DG", 1), + (construct_dg1_integral, "DG", 1), + (construct_dg2_integral, "DG", 2), pytest.param(create_dg2, "DG", 2, marks=pytest.mark.xfail(reason='Need to update TSFC in CI')), (create_cg2, "CG", 2) ]) From 87dfe977a13e16d38fe1a9ac3998f8db85366096 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 11:31:44 +0000 Subject: [PATCH 16/32] fix mistake in components --- fuse/dof.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/fuse/dof.py b/fuse/dof.py index c576ee48..97e07d50 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -190,7 +190,7 @@ def __call__(self, *args): return self.pt def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): - comps = [(i,) for v in value_shape for i in range(v)] + comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] if isinstance(self.pt, int): return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), comps if not immersed: @@ -350,6 +350,7 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): # TODO deriv dict needs implementing (currently {}) + print(value_shape) return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree, value_shape), {}, str(self)) def to_quadrature(self, arg_degree, value_shape): From ea8265b66d0ad7db60150c51c22a359ce23b65c1 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 11:32:26 +0000 Subject: [PATCH 17/32] add jacobian --- fuse/cells.py | 9 +++++++++ fuse/dof.py | 12 ++++++++---- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index b7b25833..83d426e9 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -835,6 +835,15 @@ def attachment(self, source, dst): return lambda *x: fold_reduce(attachments[0], *x) + def attachment_J_det(self, source, dst): + attachment = self.attachment(source, dst) + symbol_names = ["x", "y", "z"] + symbols = [] + for i in range(self.dim_of_node(dst)): + symbols += [sp.Symbol(symbol_names[i])] + J = sp.Matrix(attachment(*symbols)).jacobian(sp.Matrix(symbols)) + return np.sqrt(abs(float(sp.det(J.T * J)))) + def quadrature(self, degree): fiat_el = self.to_fiat() Q = create_quadrature(fiat_el, degree) diff --git a/fuse/dof.py b/fuse/dof.py index 97e07d50..3953ae33 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -375,18 +375,22 @@ def immersed(pt): if self.immersed: # need to compute jacobian from attachment. + old_pts = pts pts = np.array([self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts]) + J_det = self.cell.attachment_J_det(self.cell.id, self.cell_defined_on.id) + if not np.allclose(J_det, 1): + raise ValueError("Jacobian Determinant is not 1 did you do something wrong") # if self.pairing.orientation: # immersion = self.target_space.tabulate(wts, self.pairing.entity.orient(self.pairing.orientation))[0] # else: immersion = self.target_space.tabulate(pts, self.cell_defined_on) # Special case - force evaluation on different orientation of entity for construction of matrix transforms - if self.entity_o: - immersion = self.target_space.tabulate(wts, self.pairing.entity.orient(self.entity_o)) + # if self.entity_o: + # immersion = self.target_space.tabulate(wts, self.pairing.entity.orient(self.entity_o)) if isinstance(self.target_space, TrH1): - new_wts = wts + new_wts = wts * J_det else: - new_wts = np.outer(wts, immersion) + new_wts = np.outer(wts * J_det, immersion) else: new_wts = wts # pt dict is { pt: [(weight, component)]} From 6ad49d347b52b727399b3846da0b8c86725cd267 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 11:43:12 +0000 Subject: [PATCH 18/32] ensure vector kernel and scalar poly kernel behave the same --- fuse/dof.py | 8 +++++--- test/test_2d_examples_docs.py | 2 +- test/test_convert_to_fiat.py | 3 ++- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/fuse/dof.py b/fuse/dof.py index 3953ae33..fa02b1fb 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -190,7 +190,11 @@ def __call__(self, *args): return self.pt def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): - comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] + if len(value_shape) == 0: + comps = [[tuple()] for pt in Qpts] + else: + comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] + if isinstance(self.pt, int): return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), comps if not immersed: @@ -374,8 +378,6 @@ def immersed(pt): pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change, immersed, self.cell.dimension, value_shape) if self.immersed: - # need to compute jacobian from attachment. - old_pts = pts pts = np.array([self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts]) J_det = self.cell.attachment_J_det(self.cell.id, self.cell_defined_on.id) if not np.allclose(J_det, 1): diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index 06da1eea..1f9f1ca0 100644 --- a/test/test_2d_examples_docs.py +++ b/test/test_2d_examples_docs.py @@ -32,7 +32,7 @@ def construct_dg1(): def construct_dg0_integral(cell=None): edge = Point(1, [Point(0), Point(0)], vertex_num=2) - xs = [DOF(L2Pairing(), PolynomialKernel(1))] + xs = [DOF(L2Pairing(), VectorKernel(1))] dg0 = ElementTriple(edge, (P0, CellL2, C0), DOFGenerator(xs, S1, S1)) return dg0 diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 9bc7d3dd..9057815c 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -5,7 +5,7 @@ from firedrake import * from sympy.combinatorics import Permutation from FIAT.quadrature_schemes import create_quadrature -from test_2d_examples_docs import construct_cg1, construct_nd, construct_rt, construct_cg3, construct_dg1_integral, construct_dg2_integral +from test_2d_examples_docs import construct_cg1, construct_nd, construct_rt, construct_cg3, construct_dg0_integral, construct_dg1_integral, construct_dg2_integral from test_3d_examples_docs import construct_tet_rt, construct_tet_rt2, construct_tet_ned, construct_tet_ned_2nd_kind, construct_tet_bdm, construct_tet_ned2, construct_tet_cg4 from test_polynomial_space import flatten from element_examples import CR_n @@ -358,6 +358,7 @@ def test_entity_perms(elem_gen, cell): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg1, "CG", 1), (create_dg1, "DG", 1), + (construct_dg0_integral, "DG", 0), (construct_dg1_integral, "DG", 1), (construct_dg2_integral, "DG", 2), pytest.param(create_dg2, "DG", 2, marks=pytest.mark.xfail(reason='Need to update TSFC in CI')), From 3eaf4c72788fe58a6c293acbdd92b84046866e9c Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 12:16:45 +0000 Subject: [PATCH 19/32] bug fix --- fuse/cells.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/fuse/cells.py b/fuse/cells.py index 83d426e9..d5cbd37d 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -839,6 +839,8 @@ def attachment_J_det(self, source, dst): attachment = self.attachment(source, dst) symbol_names = ["x", "y", "z"] symbols = [] + if self.dim_of_node(dst) == 0: + return 1 for i in range(self.dim_of_node(dst)): symbols += [sp.Symbol(symbol_names[i])] J = sp.Matrix(attachment(*symbols)).jacobian(sp.Matrix(symbols)) From fdf00cc99ffaf34af66d964f1a5a762f2cc77b23 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 15:36:42 +0000 Subject: [PATCH 20/32] use ufl sobolev spaces --- fuse/spaces/element_sobolev_spaces.py | 17 ++++++++++++++++- fuse/tensor_products.py | 2 -- test/test_tensor_prod.py | 25 ++++++++++++++++--------- 3 files changed, 32 insertions(+), 12 deletions(-) diff --git a/fuse/spaces/element_sobolev_spaces.py b/fuse/spaces/element_sobolev_spaces.py index 2a3aa4fb..07c9fdd9 100644 --- a/fuse/spaces/element_sobolev_spaces.py +++ b/fuse/spaces/element_sobolev_spaces.py @@ -1,5 +1,5 @@ from functools import total_ordering - +from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2 @total_ordering class ElementSobolevSpace(object): @@ -54,6 +54,9 @@ def __init__(self, cell): def __repr__(self): return "H1" + + def to_ufl(self): + return H1 class CellHDiv(ElementSobolevSpace): @@ -64,6 +67,9 @@ def __init__(self, cell): def __repr__(self): return "HDiv" + def to_ufl(self): + return HDiv + class CellHCurl(ElementSobolevSpace): @@ -73,6 +79,9 @@ def __init__(self, cell): def __repr__(self): return "HCurl" + def to_ufl(self): + return HCurl + class CellH2(ElementSobolevSpace): @@ -82,6 +91,9 @@ def __init__(self, cell): def __repr__(self): return "H2" + def to_ufl(self): + return H2 + class CellL2(ElementSobolevSpace): @@ -90,3 +102,6 @@ def __init__(self, cell): def __repr__(self): return "L2" + + def to_ufl(self): + return L2 diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index cf16e247..4e0a486b 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -43,8 +43,6 @@ def to_ufl(self): if self.flat: return FuseElement(self, self.cell.flatten().to_ufl()) ufl_sub_elements = [e.to_ufl() for e in self.sub_elements()] - # self.setup_matrices() - # breakpoint() return TensorProductElement(*ufl_sub_elements, cell=self.cell.to_ufl()) def flatten(self): diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 8fa82f27..a90e258c 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -2,7 +2,7 @@ import numpy as np from fuse import * from firedrake import * -from test_2d_examples_docs import construct_cg1, construct_dg1 +from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg1_integral # from test_convert_to_fiat import create_cg1 @@ -32,27 +32,34 @@ def mass_solve(U): assemble(L) solve(a == L, out) assert np.allclose(out.dat.data, f.dat.data, rtol=1e-5) + return out.dat.data -@pytest.mark.parametrize("generator, code, deg", [(construct_cg1, "CG", 1), (construct_dg1, "DG", 1)]) -def test_tensor_product_ext_mesh(generator, code, deg): +@pytest.mark.parametrize("generator1, generator2, code1, code2, deg1, deg2", + [(construct_cg1, construct_cg1, "CG", "CG", 1, 1), + (construct_dg1, construct_dg1, "DG", "DG", 1, 1), + (construct_dg1, construct_cg1, "DG", "CG", 1, 1), + (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)]) +def test_ext_mesh(generator1, generator2, code1, code2, deg1, deg2): m = UnitIntervalMesh(2) mesh = ExtrudedMesh(m, 2) # manual method of creating tensor product elements - horiz_elt = FiniteElement(code, as_cell("interval"), deg) - vert_elt = FiniteElement(code, as_cell("interval"), deg) + horiz_elt = FiniteElement(code1, as_cell("interval"), deg1) + vert_elt = FiniteElement(code2, as_cell("interval"), deg2) elt = TensorProductElement(horiz_elt, vert_elt) U = FunctionSpace(mesh, elt) - mass_solve(U) + res1 = mass_solve(U) # fuseonic way of creating tensor product elements - A = generator() - B = generator() + A = generator1() + B = generator2() elem = tensor_product(A, B) U = FunctionSpace(mesh, elem.to_ufl()) - mass_solve(U) + res2 = mass_solve(U) + + assert np.allclose(res1, res2) def test_helmholtz(): From f7f27099a9c46cc79577dd560533d37fb1568b0e Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 11 Mar 2026 18:04:27 +0000 Subject: [PATCH 21/32] first stab at addition --- fuse/dof.py | 1 - fuse/spaces/polynomial_spaces.py | 4 +-- fuse/triples.py | 15 ++++++++++++ test/test_algebra.py | 42 ++++++++++++++++++++++++++++++++ 4 files changed, 58 insertions(+), 4 deletions(-) create mode 100644 test/test_algebra.py diff --git a/fuse/dof.py b/fuse/dof.py index fa02b1fb..fb04e427 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -354,7 +354,6 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): # TODO deriv dict needs implementing (currently {}) - print(value_shape) return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree, value_shape), {}, str(self)) def to_quadrature(self, arg_degree, value_shape): diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 3502a6ea..fa17c6ac 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -96,9 +96,7 @@ def __mul__(self, x): the sympy object on the right. This is due to Sympy's implementation of __mul__ not passing to this handler as it should. """ - if isinstance(x, sp.Symbol): - return ConstructedPolynomialSpace([x], [self]) - elif isinstance(x, sp.Matrix): + if isinstance(x, sp.Symbol) or isinstance(x, sp.Expr) or isinstance(x, sp.Matrix): return ConstructedPolynomialSpace([x], [self]) else: raise TypeError(f'Cannot multiply a PolySpace with {type(x)}') diff --git a/fuse/triples.py b/fuse/triples.py index 32a72208..6709eb58 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -537,6 +537,21 @@ def reverse_dof_perms(self, matrices): reversed_mats[dim][e_id] = perms_copy return reversed_mats + def __add__(self, other): + """ Construct a new element triple by combining the degrees of freedom + This implementation does not make assertions about the properties + of the resulting element. + + Elements being adding must be defined over the same cell and have the same + value shape and mapping""" + assert self.cell == other.cell + assert self.spaces[0].set_shape == other.spaces[0].set_shape + assert str(self.spaces[1]) == str(other.spaces[1]) + + spaces = (self.spaces[0] + other.spaces[0], self.spaces[1], max([self.spaces[2], other.spaces[2]])) + + return ElementTriple(self.cell, spaces, self.DOFGenerator + other.DOFGenerator) + def _to_dict(self): o_dict = {"cell": self.cell, "spaces": self.spaces, "dofs": self.DOFGenerator} return o_dict diff --git a/test/test_algebra.py b/test/test_algebra.py new file mode 100644 index 00000000..93385cab --- /dev/null +++ b/test/test_algebra.py @@ -0,0 +1,42 @@ +from fuse import * +from firedrake import * +import sympy as sp +from test_convert_to_fiat import create_cg2_tri, create_cg1, create_cr + + + +def test_bubble(): + mesh = UnitTriangleMesh() + x = SpatialCoordinate(mesh) + P2 = FiniteElement("CG", "triangle", 2) + Bubble = FiniteElement("Bubble", "triangle", 3) + P2B3 = P2 + Bubble + V = FunctionSpace(mesh, P2B3) + W = FunctionSpace(mesh, "CG", 3) + u = project(27*x[0]*x[1]*(1-x[0]-x[1]), V) + exact = Function(W) + exact.interpolate(27*x[0]*x[1]*(1-x[0]-x[1])) + # make sure that these are the same + assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14 + + +def construct_bubble(cell=None): + if cell is None: + cell = polygon(3) + x = sp.Symbol("x") + y = sp.Symbol("y") + space = PolynomialSpace(0)*(x*y*(-x-y+1)) + breakpoint() + xs = [DOF(DeltaPairing(), PointKernel((0, 0)))] + bubble = ElementTriple(cell, (space, CellL2, L2), DOFGenerator(xs, S1, S1)) + return bubble + + + +def test_temp(): + tri = polygon(3) + cg1 = create_cg1(tri) + cr1 = create_cr(tri) + cg2 = cg1 + cr1 + cg2.to_fiat() + breakpoint() \ No newline at end of file From 860f11e84854163b552774828bac15f2e0814466 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 12 Mar 2026 15:59:08 +0000 Subject: [PATCH 22/32] fix bubble - addition working --- Makefile | 6 ++++ fuse/dof.py | 4 +-- fuse/spaces/element_sobolev_spaces.py | 3 +- fuse/spaces/polynomial_spaces.py | 10 +++--- fuse/triples.py | 31 +++++++++---------- fuse/utils.py | 5 ++- test/test_algebra.py | 44 ++++++++++++--------------- 7 files changed, 51 insertions(+), 52 deletions(-) diff --git a/Makefile b/Makefile index 734a8d34..e21d433a 100644 --- a/Makefile +++ b/Makefile @@ -26,6 +26,12 @@ tests: @echo " Running all tests" @FIREDRAKE_USE_FUSE=1 python3 -m coverage run -p -m pytest -rx test +mini_tests: + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_2d_examples_docs.py + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_1d + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_orientations.py::test_surface_vec_rt + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_projection_convergence_3d\[construct_tet_ned-N1curl-1-0.8\] + coverage: @python3 -m coverage combine @python3 -m coverage report -m diff --git a/fuse/dof.py b/fuse/dof.py index fb04e427..4c6546f5 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -347,9 +347,9 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non self.pairing = self.pairing.add_entity(cell) if self.target_space is None: self.target_space = space - if self.id is None and overall_id is not None: + if overall_id is not None: self.id = overall_id - if self.sub_id is None and generator_id is not None: + if generator_id is not None: self.sub_id = generator_id def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): diff --git a/fuse/spaces/element_sobolev_spaces.py b/fuse/spaces/element_sobolev_spaces.py index 07c9fdd9..b5964d9d 100644 --- a/fuse/spaces/element_sobolev_spaces.py +++ b/fuse/spaces/element_sobolev_spaces.py @@ -1,6 +1,7 @@ from functools import total_ordering from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2 + @total_ordering class ElementSobolevSpace(object): """ @@ -54,7 +55,7 @@ def __init__(self, cell): def __repr__(self): return "H1" - + def to_ufl(self): return H1 diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index fa17c6ac..7097c5f0 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -174,9 +174,7 @@ def to_ON_polynomial_set(self, ref_el): ref_el = cell_to_simplex(ref_el) # otherwise have to work on this through tabulation - - Q = create_quadrature(ref_el, 2 * (k + 1)) - Qpts, Qwts = Q.get_points(), Q.get_weights() + weighted_sets = [] for (s, w) in zip(self.spaces, self.weights): @@ -194,6 +192,8 @@ def to_ON_polynomial_set(self, ref_el): else: vec = True w_deg = max_deg_sp_expr(w) + Q = create_quadrature(ref_el, 2 * (k + w_deg + 1)) + Qpts, Qwts = Q.get_points(), Q.get_weights() Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, shape, scale="orthonormal") # vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal") @@ -204,11 +204,11 @@ def to_ON_polynomial_set(self, ref_el): if s.set_shape or vec: scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] else: - # breakpoint() scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] scaled_at_Qpts = scaled_at_Qpts.squeeze() PkHw_coeffs = np.dot(np.multiply(scaled_at_Qpts, Qwts), Pkpw_at_Qpts.T) - # breakpoint() + if len(PkHw_coeffs.shape) == 1: + PkHw_coeffs = PkHw_coeffs.reshape(1, -1) weighted_sets.append(polynomial_set.PolynomialSet(ref_el, space.degree + w_deg, space.degree + w_deg, diff --git a/fuse/triples.py b/fuse/triples.py index 6709eb58..2d297336 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -421,7 +421,7 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): ent_dofs = entity_associations[dim][e_id][dof_gen] ent_dofs_ids = np.array([self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs], dtype=int) # (dof_gen, ent_dofs) - total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if ed.id not in total_ent_dof_ids] + total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if self.dof_id_to_fiat_id[ed.id] not in total_ent_dof_ids] # dof_idx = [total_ent_dof_ids.index(id) for id in ent_dofs_ids] dof_gen_class = ent_dofs[0].generation @@ -589,21 +589,20 @@ def num_dofs(self): return self.dof_numbers def generate(self, cell, space, id_counter): - if self.ls is None: - self.ls = [] - for l_g in self.x: - i = 0 - 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] + self.ls = [] + for l_g in self.x: + i = 0 + 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] return self.ls def make_entity_ids(self): diff --git a/fuse/utils.py b/fuse/utils.py index 6592e117..0d0ee9a9 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -67,11 +67,10 @@ def max_deg_sp_expr(sp_expr): for comp in sp_expr: # only compute degree if component is a polynomial if sp.sympify(comp).as_poly(): - degs += [sp.sympify(comp).as_poly().degree()] + degs += [sp.sympify(comp).as_poly().total_degree()] else: if sp.sympify(sp_expr).as_poly(): - degs += [sp.sympify(sp_expr).as_poly().degree()] - + degs += [sp.sympify(sp_expr).as_poly().total_degree()] return max(degs) diff --git a/test/test_algebra.py b/test/test_algebra.py index 93385cab..53049276 100644 --- a/test/test_algebra.py +++ b/test/test_algebra.py @@ -1,23 +1,8 @@ from fuse import * from firedrake import * +import numpy as np import sympy as sp -from test_convert_to_fiat import create_cg2_tri, create_cg1, create_cr - - - -def test_bubble(): - mesh = UnitTriangleMesh() - x = SpatialCoordinate(mesh) - P2 = FiniteElement("CG", "triangle", 2) - Bubble = FiniteElement("Bubble", "triangle", 3) - P2B3 = P2 + Bubble - V = FunctionSpace(mesh, P2B3) - W = FunctionSpace(mesh, "CG", 3) - u = project(27*x[0]*x[1]*(1-x[0]-x[1]), V) - exact = Function(W) - exact.interpolate(27*x[0]*x[1]*(1-x[0]-x[1])) - # make sure that these are the same - assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14 +from test_convert_to_fiat import create_cg2_tri, construct_cg3 def construct_bubble(cell=None): @@ -25,18 +10,27 @@ def construct_bubble(cell=None): cell = polygon(3) x = sp.Symbol("x") y = sp.Symbol("y") - space = PolynomialSpace(0)*(x*y*(-x-y+1)) - breakpoint() + f = (3*np.sqrt(3)/4)*(y + np.sqrt(3)/3)*(np.sqrt(3)*x + y - 2*np.sqrt(3)/3)*(-np.sqrt(3)*x + y - 2*np.sqrt(3)/3) + space = PolynomialSpace(3).restrict(0, 0)*f xs = [DOF(DeltaPairing(), PointKernel((0, 0)))] bubble = ElementTriple(cell, (space, CellL2, L2), DOFGenerator(xs, S1, S1)) return bubble +def test_bubble(): + mesh = UnitTriangleMesh() + x = SpatialCoordinate(mesh) -def test_temp(): tri = polygon(3) - cg1 = create_cg1(tri) - cr1 = create_cr(tri) - cg2 = cg1 + cr1 - cg2.to_fiat() - breakpoint() \ No newline at end of file + bub = construct_bubble(tri) + cg2 = create_cg2_tri(tri) + p2b3 = bub + cg2 + V = FunctionSpace(mesh, p2b3.to_ufl()) + W = FunctionSpace(mesh, construct_cg3().to_ufl()) + + bubble_func = 27*x[0]*x[1]*(1-x[0]-x[1]) + u = project(bubble_func, V) + exact = Function(W) + exact.interpolate(bubble_func, W) + # make sure that these are the same + assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14 From d19632043bf8ce22123dfd1257abcbda96fd9b34 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 12 Mar 2026 17:29:57 +0000 Subject: [PATCH 23/32] fix tests --- fuse/spaces/polynomial_spaces.py | 2 +- test/test_3d_examples_docs.py | 2 +- test/test_convert_to_fiat.py | 2 +- test/test_dofs.py | 2 +- test/test_perms.py | 22 +---------- test/test_tensor_prod.py | 64 ++++++++++++++++++++++++++++++-- 6 files changed, 66 insertions(+), 28 deletions(-) diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 7097c5f0..c7f8d369 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -174,7 +174,7 @@ def to_ON_polynomial_set(self, ref_el): ref_el = cell_to_simplex(ref_el) # otherwise have to work on this through tabulation - + weighted_sets = [] for (s, w) in zip(self.spaces, self.weights): diff --git a/test/test_3d_examples_docs.py b/test/test_3d_examples_docs.py index 803c3619..75729430 100644 --- a/test/test_3d_examples_docs.py +++ b/test/test_3d_examples_docs.py @@ -327,7 +327,7 @@ def test_tet_rt2(): ls = rt2.generate() # TODO make this a proper test for dof in ls: - print(dof.to_quadrature(1)) + print(dof.to_quadrature(1, value_shape=(2,))) rt2.to_fiat() diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 9057815c..ad1db363 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -358,7 +358,7 @@ def test_entity_perms(elem_gen, cell): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg1, "CG", 1), (create_dg1, "DG", 1), - (construct_dg0_integral, "DG", 0), + pytest.param(construct_dg0_integral, "DG", 0, marks=pytest.mark.xfail(reason='Passes locally, fails in CI, probably same as dg2')), (construct_dg1_integral, "DG", 1), (construct_dg2_integral, "DG", 2), pytest.param(create_dg2, "DG", 2, marks=pytest.mark.xfail(reason='Need to update TSFC in CI')), diff --git a/test/test_dofs.py b/test/test_dofs.py index ec7424f3..bf29ad45 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -220,6 +220,6 @@ def test_generate_quadrature(): print("fiat", d.pt_dict) print() for d in elem.generate(): - print("fuse", d.to_quadrature(degree)) + print("fuse", d.to_quadrature(degree, value_shape=(2,))) elem.to_fiat() diff --git a/test/test_perms.py b/test/test_perms.py index f00aed6e..82a3cb90 100644 --- a/test/test_perms.py +++ b/test/test_perms.py @@ -1,9 +1,8 @@ from fuse import * -from test_convert_to_fiat import create_cg1, create_dg1, create_cg2 +from test_convert_to_fiat import create_cg1, create_dg1, create_cg2, construct_nd from test_2d_examples_docs import construct_cg3 import pytest import numpy as np -import sympy as sp vert = Point(0) edge = Point(1, [Point(0), Point(0)], vertex_num=2) @@ -48,24 +47,7 @@ def test_basic_perms(cell): @pytest.mark.parametrize("cell", [tri]) def test_nd_perms(cell): - deg = 1 - edge = cell.edges(get_class=True)[0] - x = sp.Symbol("x") - y = sp.Symbol("y") - - xs = [DOF(L2Pairing(), PolynomialKernel(edge.basis_vectors()[0]))] - dofs = DOFGenerator(xs, S1, S2) - int_ned = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - - xs = [immerse(cell, int_ned, TrHCurl)] - tri_dofs = DOFGenerator(xs, C3, S3) - - M = sp.Matrix([[y, -x]]) - vec_Pk = PolynomialSpace(deg - 1, set_shape=True) - Pk = PolynomialSpace(deg - 1) - nd = vec_Pk + (Pk.restrict(deg - 2, deg - 1))*M - - ned = ElementTriple(cell, (nd, CellHCurl, C0), [tri_dofs]) + ned = construct_nd(cell) ned.to_fiat() for i, mat in ned.matrices[2][0].items(): print(i) diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index a90e258c..15e15a71 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -36,10 +36,10 @@ def mass_solve(U): @pytest.mark.parametrize("generator1, generator2, code1, code2, deg1, deg2", - [(construct_cg1, construct_cg1, "CG", "CG", 1, 1), - (construct_dg1, construct_dg1, "DG", "DG", 1, 1), - (construct_dg1, construct_cg1, "DG", "CG", 1, 1), - (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)]) + [(construct_cg1, construct_cg1, "CG", "CG", 1, 1), + (construct_dg1, construct_dg1, "DG", "DG", 1, 1), + (construct_dg1, construct_cg1, "DG", "CG", 1, 1), + (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)]) def test_ext_mesh(generator1, generator2, code1, code2, deg1, deg2): m = UnitIntervalMesh(2) mesh = ExtrudedMesh(m, 2) @@ -137,3 +137,59 @@ def test_flattening(A, B, res): else: cell = tensor_cell.flatten() cell.construct_fuse_rep() + + +# def test_trace_galerkin_projection(): +# mesh = UnitSquareMesh(10, 10, quadrilateral=True) + +# x, y = SpatialCoordinate(mesh) +# A = construct_cg1() +# B = construct_dg0_integral() +# elem = tensor_product(A, B) +# elem = elem.flatten() + +# # Define the Trace Space +# T = FunctionSpace(mesh, elem.to_ufl()) + +# # Define trial and test functions +# lambdar = TrialFunction(T) +# gammar = TestFunction(T) + +# # Define right hand side function + +# V = FunctionSpace(mesh, "CG", 1) +# f = Function(V) +# f.interpolate(cos(x*pi*2)*cos(y*pi*2)) + +# # Construct bilinear form +# a = inner(lambdar, gammar) * ds + inner(lambdar('+'), gammar('+')) * dS + +# # Construct linear form +# l = inner(f, gammar) * ds + inner(f('+'), gammar('+')) * dS + +# # Compute the solution +# t = Function(T) +# solve(a == l, t, solver_parameters={'ksp_rtol': 1e-14}) + +# # Compute error in trace norm +# trace_error = sqrt(assemble(FacetArea(mesh)*inner((t - f)('+'), (t - f)('+')) * dS)) + +# assert trace_error < 1e-13 + +# def test_hdiv(): +# np.set_printoptions(linewidth=90, precision=4, suppress=True) +# m = UnitIntervalMesh(2) +# mesh = ExtrudedMesh(m, 2) +# CG_1 = FiniteElement("CG", interval, 1) +# DG_0 = FiniteElement("DG", interval, 0) +# P1P0 = TensorProductElement(CG_1, DG_0) +# RT_horiz = HDivElement(P1P0) +# P0P1 = TensorProductElement(DG_0, CG_1) +# RT_vert = HDivElement(P0P1) +# elt = RT_horiz + RT_vert +# V = FunctionSpace(mesh, elt) +# tabulation = V.finat_element.fiat_equivalent.tabulate(0, [(0, 0), (1, 0)]) +# for ent, arr in tabulation.items(): +# print(ent) +# for comp in arr: +# print(comp[0], comp[1]) From ec273f63a639b51305f09b6ded517b9e982c8dea Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 19 Mar 2026 16:48:53 +0000 Subject: [PATCH 24/32] working towards symmetric tensor products --- fuse/cells.py | 9 ++++ fuse/tensor_products.py | 58 +++++++++++++++++++--- fuse/triples.py | 32 +++++++----- test/test_convert_to_fiat.py | 4 +- test/test_tensor_prod.py | 96 ++++++++++++++++++++---------------- 5 files changed, 137 insertions(+), 62 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index d5cbd37d..4fc292bf 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -980,6 +980,15 @@ def __init__(self, A, B): def ordered_vertices(self): return self.A.ordered_vertices() + self.B.ordered_vertices() + def get_starter_ids(self): + # this doesn't actually make sense - remove when confirmed all changes to eliminate min ids from triple is done + a_starts = self.A.get_starter_ids() + b_starts = self.B.get_starter_ids() + ids = [] + for a, b in zip(a_starts, b_starts): + ids += [max(a, b)] + return ids + def get_spatial_dimension(self): return self.dimension diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index 4e0a486b..21452353 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -22,6 +22,10 @@ def __init__(self, A, B, flat=False): self.cell = TensorProductPoint(A.cell, B.cell) self.flat = flat self.apply_matrices = False + if self.flat: + self.unflat_cell = self.cell + self.cell = self.cell.flatten() + self.setup_matrices() def sub_elements(self): return [self.A, self.B] @@ -30,18 +34,58 @@ def __repr__(self): return "TensorProd(%s, %s)" % (repr(self.A), repr(self.B)) def setup_matrices(self): - oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(self.A.generate() + self.B.generate()) + if self.A.dimension > 1 or self.B.dimension > 1: + raise NotImplementedError("Combining of matrices not implemented in 3D") + self.A.to_ufl() + self.B.to_ufl() + oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(self.generate()) + unflat_top = self.unflat_cell.to_fiat().get_topology() + for dim in unflat_top.keys(): + a_ents = self.A.cell.get_topology()[dim[0]].keys() + b_ents = self.A.cell.get_topology()[dim[1]].keys() + ents = [(a, b) for a in a_ents for b in b_ents] + for e, (a, b) in enumerate(ents): + ent_dofs = self.entity_dofs[dim][(a, b)] + if len(ent_dofs) > 1: + sub_mat = oriented_mats_by_entity[sum(dim)][e] + a_mat = self.A.matrices[dim[0]][a] + b_mat = self.B.matrices[dim[1]][b] + # need to make groups for tensor product cell that are + # different for flat or not + breakpoint() + from collections import defaultdict + from FIAT.reference_element import tuple_sum + breakpoint() - for dim in range(self.cell.dimension): - for dimA in range(self.A.cell.dimension): - pass - for dimB in range(self.B.cell_dimension): - pass return super().setup_matrices() + def generate(self): + a_dofs = self.A.generate() + b_dofs = self.B.generate() + a_ent_assocs, _, _ = self.A._entity_associations(a_dofs, overall=False) + b_ent_assocs, _, _ = self.B._entity_associations(b_dofs, overall=False) + unflat_top = self.unflat_cell.to_fiat().get_topology() + self.entity_dofs = {} + dofs = [] + counter = 0 + for dim in unflat_top.keys(): + ents_A = a_ent_assocs[dim[0]].keys() + ents_B = b_ent_assocs[dim[1]].keys() + self.entity_dofs[dim] = {(a_e, b_e): tuple() for a_e in ents_A for b_e in ents_B} + for a_e, b_e in self.entity_dofs[dim].keys(): + a_dofs = [d for dofs in a_ent_assocs[dim[0]][a_e].values() for d in dofs] + b_dofs = [d for dofs in b_ent_assocs[dim[1]][b_e].values() for d in dofs] + new_dofs = [(a, b) for a in a_dofs for b in b_dofs] + dofs += new_dofs + self.entity_dofs[dim][(a_e, b_e)] = [i + counter for i in range(len(new_dofs))] + counter += len(new_dofs) + self.dofs = dofs + return dofs + + def to_ufl(self): if self.flat: - return FuseElement(self, self.cell.flatten().to_ufl()) + return FuseElement(self, self.cell.to_ufl()) ufl_sub_elements = [e.to_ufl() for e in self.sub_elements()] return TensorProductElement(*ufl_sub_elements, cell=self.cell.to_ufl()) diff --git a/fuse/triples.py b/fuse/triples.py index 2d297336..6c220383 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -343,9 +343,8 @@ def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): res_dict[dim][e_id][val] = np.matmul(transformed_basis, original_V.T) return res_dict - def _entity_associations(self, dofs): - min_ids = self.cell.get_starter_ids() - entity_associations = {dim: {e.id - min_ids[dim]: {} for e in self.cell.d_entities(dim)} + def _entity_associations(self, dofs, overall=True): + entity_associations = {dim: {i: {} for i, e in enumerate(self.cell.d_entities(dim))} for dim in range(self.cell.dim() + 1)} cell_dim = self.cell.dim() cell_dict = entity_associations[cell_dim][0] @@ -355,8 +354,13 @@ def _entity_associations(self, dofs): # construct mapping of entities to the dof generators and the dofs they generate for d in dofs: sub_dim = d.cell_defined_on.dim() - sub_dict = entity_associations[sub_dim][d.cell_defined_on.id - min_ids[sub_dim]] - for dim in set([sub_dim, cell_dim]): + cell_defined_on_id = self.cell.d_entities_ids(sub_dim).index(d.cell_defined_on.id) + sub_dict = entity_associations[sub_dim][cell_defined_on_id] + if overall: + dims = set([sub_dim, cell_dim]) + else: + dims = [sub_dim] + for dim in dims: dof_gen = str(d.generation[dim]) if not len(d.generation[dim].g2.members()) == 1: @@ -369,7 +373,6 @@ def _entity_associations(self, dofs): sub_dict[dof_gen] += [d] elif dim < cell_dim or not d.immersed: sub_dict[dof_gen] = [d] - if dof_gen in cell_dict.keys() and dim == cell_dim and d.immersed: cell_dict[dof_gen] += [d] elif dim == cell_dim and d.immersed: @@ -385,8 +388,10 @@ def _initialise_entity_dicts(self, dofs): oriented_mats_by_entity[dim] = {} flat_by_entity[dim] = {} ents = self.cell.d_entities(dim) - for e in ents: - e_id = e.id - min_ids[dim] + for e_id, e in enumerate(ents): + old_e_id = e.id - min_ids[dim] + if e_id != old_e_id: + raise ValueError("e id problems") members = e.group.members() oriented_mats_by_entity[dim][e_id] = {} flat_by_entity[dim][e_id] = {} @@ -411,8 +416,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): # dof mapping according to the generation for dim in range(self.cell.dim() + 1): ents = self.cell.d_entities(dim) - for e in ents: - e_id = e.id - min_ids[dim] + for e_id, e in enumerate(ents): + old_e_id = e.id - min_ids[dim] + if e_id != old_e_id: + raise ValueError("e id problems") members = e.group.members() for g in members: val = g.numeric_rep() @@ -478,7 +485,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): g_sub_mat = perm_list_to_matrix(identity, [sub_e for sub_e, _ in permuted_ents]) for sub_e, sub_g in permuted_ents: sub_e = self.cell.get_node(sub_e) - sub_e_id = sub_e.id - min_ids[sub_e.dim()] + old_sub_e_id = sub_e.id - min_ids[sub_e.dim()] + sub_e_id = self.cell.d_entities(sub_e.dim(), get_class=False).index(sub_e.id) + if sub_e_id != old_sub_e_id: + raise ValueError("sub e id problems") sub_ent_ids = [] for (k, v) in entity_associations[immersed_dim][sub_e_id].items(): sub_ent_ids += [self.dof_id_to_fiat_id[e.id] for e in v] diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index ad1db363..01f5b692 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -140,7 +140,9 @@ def create_cg1_flipped(cell): return cg -def create_cg2(cell): +def create_cg2(cell=None): + if cell == None: + cell = line() deg = 2 if cell.dim() > 1: raise NotImplementedError("This method is for cg2 on edges, please use create_cg2_tri for triangles") diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 15e15a71..736732dc 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -2,7 +2,8 @@ import numpy as np from fuse import * from firedrake import * -from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg1_integral +from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg0_integral, construct_dg1_integral +from test_convert_to_fiat import create_cg2 # from test_convert_to_fiat import create_cg1 @@ -96,34 +97,26 @@ def test_on_quad_mesh(): U = FunctionSpace(m, "CG", 1) mass_solve(U) - -def test_quad_mesh_helmholtz(): +@pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(construct_cg1, "CG", 1, 1.8), + (create_cg2, "CG", 2, 3.8)]) +def test_quad_mesh_helmholtz(elem_gen, elem_code, deg, conv_rate): quadrilateral = True vals = range(3, 6) res_fuse = [] - res_fire = [] for r in vals: mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) - A = construct_cg1() - B = construct_cg1() + A = elem_gen() + B = elem_gen() elem = tensor_product(A, B).flatten() U = FunctionSpace(mesh, elem.to_ufl()) res_fuse += [helmholtz_solve(mesh, U)] - U = FunctionSpace(mesh, "CG", 1) - res_fire += [helmholtz_solve(mesh, U)] - print("l2 error norms:", res_fuse) + print("Fuse l2 error norms:", res_fuse) res = np.array(res_fuse) conv = np.log2(res[:-1] / res[1:]) - print("convergence order:", conv) - assert (np.array(conv) > 1.8).all() - - print("l2 error norms:", res_fire) - res = np.array(res_fire) - conv = np.log2(res[:-1] / res[1:]) - print("convergence order:", conv) - assert (np.array(conv) > 1.8).all() + print("Fuse convergence order:", conv) + assert (np.array(conv) > conv_rate).all() @pytest.mark.parametrize(["A", "B", "res"], [(Point(0), line(), False), @@ -138,43 +131,60 @@ def test_flattening(A, B, res): cell = tensor_cell.flatten() cell.construct_fuse_rep() +def test_cg1_dg0(): + A = construct_cg1() + B = construct_dg1_integral() + non_sym = tensor_product(A, B).flatten() + non_sym2 = tensor_product(B, A).flatten() + # from finat.element_factory import convert + # non_sym, _ = convert(non_sym.to_ufl(), shift_axes=0) + # non_sym2, _ = convert(non_sym2.to_ufl(), shift_axes=0) + # from FIAT.reference_element import flatten_entities + # print() + # print(non_sym.entity_dofs()) + # print(non_sym2.entity_dofs()) + # print(flatten_entities(non_sym.entity_dofs())) + # print(flatten_entities(non_sym2.entity_dofs())) + breakpoint() + + +def test_trace_galerkin_projection(): + mesh = UnitSquareMesh(10, 10, quadrilateral=True) -# def test_trace_galerkin_projection(): -# mesh = UnitSquareMesh(10, 10, quadrilateral=True) + x, y = SpatialCoordinate(mesh) + A = construct_cg1() + B = construct_dg1_integral() + elem = tensor_product(A, B).flatten() + elem2 = tensor_product(B, A).flatten() -# x, y = SpatialCoordinate(mesh) -# A = construct_cg1() -# B = construct_dg0_integral() -# elem = tensor_product(A, B) -# elem = elem.flatten() -# # Define the Trace Space -# T = FunctionSpace(mesh, elem.to_ufl()) + # Define the Trace Space + T = FunctionSpace(mesh, elem.to_ufl() + elem2.to_ufl()) -# # Define trial and test functions -# lambdar = TrialFunction(T) -# gammar = TestFunction(T) + # Define trial and test functions + lambdar = TrialFunction(T) + gammar = TestFunction(T) -# # Define right hand side function + # Define right hand side function -# V = FunctionSpace(mesh, "CG", 1) -# f = Function(V) -# f.interpolate(cos(x*pi*2)*cos(y*pi*2)) + V = FunctionSpace(mesh, "CG", 1) + f = Function(V) + f.interpolate(cos(x*pi*2)*cos(y*pi*2)) -# # Construct bilinear form -# a = inner(lambdar, gammar) * ds + inner(lambdar('+'), gammar('+')) * dS + # Construct bilinear form + a = inner(lambdar, gammar) * ds + inner(lambdar('+'), gammar('+')) * dS -# # Construct linear form -# l = inner(f, gammar) * ds + inner(f('+'), gammar('+')) * dS + # Construct linear form + l = inner(f, gammar) * ds + inner(f('+'), gammar('+')) * dS -# # Compute the solution -# t = Function(T) -# solve(a == l, t, solver_parameters={'ksp_rtol': 1e-14}) + # Compute the solution + t = Function(T) + solve(a == l, t, solver_parameters={'ksp_rtol': 1e-14}) -# # Compute error in trace norm -# trace_error = sqrt(assemble(FacetArea(mesh)*inner((t - f)('+'), (t - f)('+')) * dS)) + # Compute error in trace norm + trace_error = sqrt(assemble(FacetArea(mesh)*inner((t - f)('+'), (t - f)('+')) * dS)) -# assert trace_error < 1e-13 + assert trace_error < 1e-13 # def test_hdiv(): # np.set_printoptions(linewidth=90, precision=4, suppress=True) From f3627a0698599127bedf40db708c7566d3289fdd Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 20 Mar 2026 16:09:50 +0000 Subject: [PATCH 25/32] add groups to non flattened tensor product cell --- fuse/cells.py | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 4fc292bf..24a49ba8 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -976,10 +976,25 @@ def __init__(self, A, B): self.B = B self.dimension = self.A.dimension + self.B.dimension self.flat = False + self.fiat_elem = None + self.group = self.compute_cell_group().add_cell(self) def ordered_vertices(self): return self.A.ordered_vertices() + self.B.ordered_vertices() + def compute_cell_group(self): + """ + Systematically work out the symmetry group of the tensor product cell. + """ + verts = self.vertices() + group = [(g_a, g_b) for g_a in self.A.group.members() for g_b in self.A.group.members()] + perms = [] + for g_a, g_b in group: + new_verts = [(v_a, v_b) for v_a in g_a.permute(self.A.vertices()) for v_b in g_b.permute(self.B.vertices())] + perm = [verts.index(v) for v in new_verts] + perms += [fuse_groups.Permutation(perm)] + return fuse_groups.PermutationSetRepresentation(perms) + def get_starter_ids(self): # this doesn't actually make sense - remove when confirmed all changes to eliminate min ids from triple is done a_starts = self.A.get_starter_ids() @@ -993,8 +1008,7 @@ def get_spatial_dimension(self): return self.dimension def get_sub_entities(self): - self.A.get_sub_entities() - self.B.get_sub_entities() + return self.to_fiat().sub_entities def dim(self): return (self.A.dimension, self.B.dimension) @@ -1004,18 +1018,19 @@ def d_entities(self, d, get_class=True): def vertices(self, get_class=True, return_coords=False): # TODO maybe refactor with get_node - verts = self.d_entities(0, get_class) if return_coords: a_verts = self.A.vertices(return_coords=return_coords) b_verts = self.B.vertices(return_coords=return_coords) return [a + b for a in a_verts for b in b_verts] - return verts + return [(a, b) for a in self.A.vertices() for b in self.B.vertices()] def to_ufl(self, name=None): return TensorProductCell(self.A.to_ufl(), self.B.to_ufl()) def to_fiat(self, name=None): - return CellComplexToFiatTensorProduct(self, name) + if self.fiat_elem is None: + self.fiat_elem = CellComplexToFiatTensorProduct(self, name) + return self.fiat_elem def flatten(self): assert self.A.equivalent(self.B) From 8f3719a56c0d1970d1e874209e593de21bb6de0b Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 23 Mar 2026 12:16:16 +0000 Subject: [PATCH 26/32] flattening or not --- fuse/tensor_products.py | 21 ++++++++++++++------- fuse/triples.py | 6 ++---- test/test_tensor_prod.py | 4 ++-- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index 21452353..028dfc93 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -22,10 +22,10 @@ def __init__(self, A, B, flat=False): self.cell = TensorProductPoint(A.cell, B.cell) self.flat = flat self.apply_matrices = False + self.setup_matrices() if self.flat: self.unflat_cell = self.cell self.cell = self.cell.flatten() - self.setup_matrices() def sub_elements(self): return [self.A, self.B] @@ -34,13 +34,18 @@ def __repr__(self): return "TensorProd(%s, %s)" % (repr(self.A), repr(self.B)) def setup_matrices(self): - if self.A.dimension > 1 or self.B.dimension > 1: + if self.A.cell.dimension > 1 or self.B.cell.dimension > 1: raise NotImplementedError("Combining of matrices not implemented in 3D") + if self.cell.flat: + raise NotImplementedError("Matrices for flattened cells not yet implemented") self.A.to_ufl() self.B.to_ufl() oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(self.generate()) - unflat_top = self.unflat_cell.to_fiat().get_topology() - for dim in unflat_top.keys(): + if self.flat: + top = self.unflat_cell.to_fiat().get_topology() + else: + top = self.cell.to_fiat().get_topology() + for dim in top.keys(): a_ents = self.A.cell.get_topology()[dim[0]].keys() b_ents = self.A.cell.get_topology()[dim[1]].keys() ents = [(a, b) for a in a_ents for b in b_ents] @@ -64,11 +69,14 @@ def generate(self): b_dofs = self.B.generate() a_ent_assocs, _, _ = self.A._entity_associations(a_dofs, overall=False) b_ent_assocs, _, _ = self.B._entity_associations(b_dofs, overall=False) - unflat_top = self.unflat_cell.to_fiat().get_topology() + if self.flat: + top = self.unflat_cell.to_fiat().get_topology() + else: + top = self.cell.to_fiat().get_topology() self.entity_dofs = {} dofs = [] counter = 0 - for dim in unflat_top.keys(): + for dim in top.keys(): ents_A = a_ent_assocs[dim[0]].keys() ents_B = b_ent_assocs[dim[1]].keys() self.entity_dofs[dim] = {(a_e, b_e): tuple() for a_e in ents_A for b_e in ents_B} @@ -82,7 +90,6 @@ def generate(self): self.dofs = dofs return dofs - def to_ufl(self): if self.flat: return FuseElement(self, self.cell.to_ufl()) diff --git a/fuse/triples.py b/fuse/triples.py index 6c220383..b7b606ef 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -384,14 +384,12 @@ def _initialise_entity_dicts(self, dofs): dof_id_mat = np.eye(len(dofs)) oriented_mats_by_entity = {} flat_by_entity = {} - for dim in range(self.cell.dim() + 1): + for dim in range(self.cell.dimension + 1): oriented_mats_by_entity[dim] = {} flat_by_entity[dim] = {} ents = self.cell.d_entities(dim) for e_id, e in enumerate(ents): - old_e_id = e.id - min_ids[dim] - if e_id != old_e_id: - raise ValueError("e id problems") + # old_e_id = e.id - min_ids[dim] members = e.group.members() oriented_mats_by_entity[dim][e_id] = {} flat_by_entity[dim][e_id] = {} diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 736732dc..260f84cd 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -134,8 +134,8 @@ def test_flattening(A, B, res): def test_cg1_dg0(): A = construct_cg1() B = construct_dg1_integral() - non_sym = tensor_product(A, B).flatten() - non_sym2 = tensor_product(B, A).flatten() + non_sym = tensor_product(A, B) + # non_sym2 = tensor_product(B, A).flatten() # from finat.element_factory import convert # non_sym, _ = convert(non_sym.to_ufl(), shift_axes=0) # non_sym2, _ = convert(non_sym2.to_ufl(), shift_axes=0) From 6e4f40fbe8564a443bb04886bb34b71d770993c7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 20 Apr 2026 20:51:34 +0100 Subject: [PATCH 27/32] work on combining matrices --- fuse/cells.py | 33 +++++++++++++++++++--- fuse/tensor_products.py | 26 +++++++++++------ fuse/triples.py | 16 ++++++++--- test/test_convert_to_fiat.py | 2 +- test/test_tensor_prod.py | 55 +++++++++++++++++++++++------------- 5 files changed, 95 insertions(+), 37 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 24a49ba8..fda3f194 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -972,28 +972,50 @@ def _from_dict(o_dict): class TensorProductPoint(): def __init__(self, A, B): + self.id = next(self.id_iter) self.A = A self.B = B self.dimension = self.A.dimension + self.B.dimension self.flat = False self.fiat_elem = None - self.group = self.compute_cell_group().add_cell(self) + self.group = self.compute_cell_group() def ordered_vertices(self): return self.A.ordered_vertices() + self.B.ordered_vertices() + def component_orientations(self): + from fuse.utils import orientation_value + self.component_os_to_os = {} + for dim in self.to_fiat().get_topology(): + self.component_os_to_os[dim] = {} + a_ent = self.A.d_entities(dim[0])[0] + b_ent = self.B.d_entities(dim[1])[0] + verts = [(v_a, v_b) for v_a in a_ent.vertices() for v_b in b_ent.vertices()] + ident = [i for i in range(len(verts))] + group = [(g_a, g_b) for g_a in a_ent.group.members() for g_b in b_ent.group.members()] + for g_a, g_b in group: + new_verts = [(v_a, v_b) for v_a in g_a.permute(a_ent.vertices()) for v_b in g_b.permute(b_ent.vertices())] + perm = [verts.index(v) for v in new_verts] + o_val = orientation_value(ident, perm) + if sum(dim) == self.dimension and self.group.group_rep_numbering is not None: + o_val = self.group.group_rep_numbering[o_val] + self.component_os_to_os[dim][(g_a.numeric_rep(), g_b.numeric_rep())] = o_val + return self.component_os_to_os + def compute_cell_group(self): """ Systematically work out the symmetry group of the tensor product cell. """ verts = self.vertices() - group = [(g_a, g_b) for g_a in self.A.group.members() for g_b in self.A.group.members()] + group = [(g_a, g_b) for g_a in self.A.group.members() for g_b in self.B.group.members()] perms = [] for g_a, g_b in group: new_verts = [(v_a, v_b) for v_a in g_a.permute(self.A.vertices()) for v_b in g_b.permute(self.B.vertices())] perm = [verts.index(v) for v in new_verts] perms += [fuse_groups.Permutation(perm)] - return fuse_groups.PermutationSetRepresentation(perms) + + grp = fuse_groups.PermutationSetRepresentation(perms).add_cell(self) + return grp def get_starter_ids(self): # this doesn't actually make sense - remove when confirmed all changes to eliminate min ids from triple is done @@ -1014,6 +1036,9 @@ def dim(self): return (self.A.dimension, self.B.dimension) def d_entities(self, d, get_class=True): + if isinstance(d, tuple): + return [TensorProductPoint(e_a, e_b) for e_a in self.A.d_entities(d[0], get_class) for e_b in self.B.d_entities(d[1], get_class)] + raise NotImplementedError("not sure this is right") return self.A.d_entities(d, get_class) + self.B.d_entities(d, get_class) def vertices(self, get_class=True, return_coords=False): @@ -1029,7 +1054,7 @@ def to_ufl(self, name=None): def to_fiat(self, name=None): if self.fiat_elem is None: - self.fiat_elem = CellComplexToFiatTensorProduct(self, name) + self.fiat_elem = CellComplexToFiatTensorProduct(self, name) return self.fiat_elem def flatten(self): diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index 028dfc93..8891ce48 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -1,5 +1,6 @@ from fuse.triples import ElementTriple from fuse.cells import TensorProductPoint +import numpy as np from finat.ufl import TensorProductElement, FuseElement @@ -47,22 +48,29 @@ def setup_matrices(self): top = self.cell.to_fiat().get_topology() for dim in top.keys(): a_ents = self.A.cell.get_topology()[dim[0]].keys() - b_ents = self.A.cell.get_topology()[dim[1]].keys() + b_ents = self.B.cell.get_topology()[dim[1]].keys() ents = [(a, b) for a in a_ents for b in b_ents] + comp_os = self.cell.component_orientations() for e, (a, b) in enumerate(ents): ent_dofs = self.entity_dofs[dim][(a, b)] if len(ent_dofs) > 1: - sub_mat = oriented_mats_by_entity[sum(dim)][e] + sub_mat = oriented_mats_by_entity[dim][e] a_mat = self.A.matrices[dim[0]][a] + a_ent_ids = self.A.entity_ids[dim[0]][a] b_mat = self.B.matrices[dim[1]][b] - # need to make groups for tensor product cell that are - # different for flat or not - breakpoint() - from collections import defaultdict - from FIAT.reference_element import tuple_sum + b_ent_ids = self.B.entity_ids[dim[1]][b] - breakpoint() - return super().setup_matrices() + os = [(o_a, o_b) for o_a in a_mat.keys() for o_b in b_mat.keys()] + for o in os: + a_sub_mat = a_mat[o[0]][np.ix_(a_ent_ids, a_ent_ids)] + b_sub_mat = b_mat[o[1]][np.ix_(b_ent_ids, b_ent_ids)] + combined_sub_mat = np.kron(a_sub_mat, b_sub_mat) + new_o = comp_os[dim][o] + sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)] = np.matmul(sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)], combined_sub_mat) + # from collections import defaultdict + # from FIAT.reference_element import tuple_sum + self.matrices = oriented_mats_by_entity + self.reversed_matrices = self.reverse_dof_perms(self.matrices) def generate(self): a_dofs = self.A.generate() diff --git a/fuse/triples.py b/fuse/triples.py index b7b606ef..4a68e544 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -380,13 +380,18 @@ def _entity_associations(self, dofs, overall=True): return entity_associations, pure_perm, sub_pure_perm def _initialise_entity_dicts(self, dofs): - min_ids = self.cell.get_starter_ids() + # min_ids = self.cell.get_starter_ids() dof_id_mat = np.eye(len(dofs)) oriented_mats_by_entity = {} flat_by_entity = {} - for dim in range(self.cell.dimension + 1): + if isinstance(self.cell, TensorProductPoint): + dims = [(a_d, b_d) for a_d in range(self.cell.A.dimension + 1) for b_d in range(self.cell.B.dimension + 1)] + else: + dims = [i for i in range(self.cell.dimension + 1)] + for dim in dims: oriented_mats_by_entity[dim] = {} flat_by_entity[dim] = {} + ents = self.cell.d_entities(dim) for e_id, e in enumerate(ents): # old_e_id = e.id - min_ids[dim] @@ -531,13 +536,16 @@ def orient_mat_perms(self): num_ents += len(ents) def reverse_dof_perms(self, matrices): - min_ids = self.cell.get_starter_ids() + # min_ids = self.cell.get_starter_ids() reversed_mats = {} for dim in matrices.keys(): reversed_mats[dim] = {} ents = self.cell.d_entities(dim) for e in ents: - e_id = e.id - min_ids[dim] + # old_e_id = e.id - min_ids[dim] + e_id = self.cell.d_entities(e.dim(), get_class=False).index(e.id) + if isinstance(dim, tuple): + breakpoint() perms_copy = matrices[dim][e_id].copy() members = e.group.members() for m in members: diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 01f5b692..1a655762 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -141,7 +141,7 @@ def create_cg1_flipped(cell): def create_cg2(cell=None): - if cell == None: + if cell is None: cell = line() deg = 2 if cell.dim() > 1: diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 260f84cd..7fbb158b 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -2,7 +2,7 @@ import numpy as np from fuse import * from firedrake import * -from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg0_integral, construct_dg1_integral +from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg1_integral from test_convert_to_fiat import create_cg2 # from test_convert_to_fiat import create_cg1 @@ -97,6 +97,7 @@ def test_on_quad_mesh(): U = FunctionSpace(m, "CG", 1) mass_solve(U) + @pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(construct_cg1, "CG", 1, 1.8), (create_cg2, "CG", 2, 3.8)]) def test_quad_mesh_helmholtz(elem_gen, elem_code, deg, conv_rate): @@ -131,6 +132,7 @@ def test_flattening(A, B, res): cell = tensor_cell.flatten() cell.construct_fuse_rep() + def test_cg1_dg0(): A = construct_cg1() B = construct_dg1_integral() @@ -145,6 +147,7 @@ def test_cg1_dg0(): # print(non_sym2.entity_dofs()) # print(flatten_entities(non_sym.entity_dofs())) # print(flatten_entities(non_sym2.entity_dofs())) + print(non_sym) breakpoint() @@ -157,7 +160,6 @@ def test_trace_galerkin_projection(): elem = tensor_product(A, B).flatten() elem2 = tensor_product(B, A).flatten() - # Define the Trace Space T = FunctionSpace(mesh, elem.to_ufl() + elem2.to_ufl()) @@ -186,20 +188,35 @@ def test_trace_galerkin_projection(): assert trace_error < 1e-13 -# def test_hdiv(): -# np.set_printoptions(linewidth=90, precision=4, suppress=True) -# m = UnitIntervalMesh(2) -# mesh = ExtrudedMesh(m, 2) -# CG_1 = FiniteElement("CG", interval, 1) -# DG_0 = FiniteElement("DG", interval, 0) -# P1P0 = TensorProductElement(CG_1, DG_0) -# RT_horiz = HDivElement(P1P0) -# P0P1 = TensorProductElement(DG_0, CG_1) -# RT_vert = HDivElement(P0P1) -# elt = RT_horiz + RT_vert -# V = FunctionSpace(mesh, elt) -# tabulation = V.finat_element.fiat_equivalent.tabulate(0, [(0, 0), (1, 0)]) -# for ent, arr in tabulation.items(): -# print(ent) -# for comp in arr: -# print(comp[0], comp[1]) + +def test_hdiv(): + np.set_printoptions(linewidth=90, precision=4, suppress=True) + m = UnitIntervalMesh(2) + mesh = ExtrudedMesh(m, 2) + CG_1 = FiniteElement("CG", "interval", 1) + DG_0 = FiniteElement("DG", "interval", 0) + P1P0 = TensorProductElement(CG_1, DG_0) + RT_horiz = HDivElement(P1P0) + P0P1 = TensorProductElement(DG_0, CG_1) + RT_vert = HDivElement(P0P1) + elt = RT_horiz + RT_vert + # mesh = UnitSquareMesh(1, 1, quadrilateral=True) + V = FunctionSpace(mesh, elt) + u = TrialFunction(V) + v = TestFunction(V) + f = Function(V) + x, y = SpatialCoordinate(mesh) + f_vec = as_vector(((1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2), (1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2))) + f = project(f_vec, V) + a = (inner(grad(u), grad(v)) + inner(u, v)) * dx + L = inner(f, v) * dx + u = Function(V) + solve(a == L, u) + breakpoint() + # f.interpolate(cos(x*pi*2)*cos(y*pi*2)) + # V.finat_element.basis_evaluation(1, [(0, 0)]) + # tabulation = V.finat_element.fiat_equivalent.tabulate(0, [(0, 0), (1, 0)]) + # for ent, arr in tabulation.items(): + # print(ent) + # for comp in arr: + # print(comp[0], comp[1]) From f91d745c4f5f09a8f6a86229647cc857a527e0f3 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 21 Apr 2026 19:36:46 +0100 Subject: [PATCH 28/32] fix choice of form degree for fuse elements, work on making hdiv element take fuse elements --- fuse/cells.py | 14 +++++++++++++- fuse/dof.py | 2 +- fuse/triples.py | 19 ++++++++++++++++--- test/test_2d_examples_docs.py | 12 +++++++++++- test/test_tensor_prod.py | 19 +++++++++++++++---- 5 files changed, 56 insertions(+), 10 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index fda3f194..b7a9d1d4 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -970,6 +970,7 @@ def _from_dict(o_dict): class TensorProductPoint(): + id_iter = itertools.count() def __init__(self, A, B): self.id = next(self.id_iter) @@ -979,6 +980,15 @@ def __init__(self, A, B): self.flat = False self.fiat_elem = None self.group = self.compute_cell_group() + self.entities = {} + for d in [(a_d, b_d) for a_d in range(self.A.dimension + 1) for b_d in range(self.B.dimension + 1)]: + if d == (self.A.dimension, self.B.dimension): + self.entities[d] = [self] + else: + self.entities[d] = [TensorProductPoint(e_a, e_b) for e_a in self.A.d_entities(d[0], True) for e_b in self.B.d_entities(d[1], True)] + + + def ordered_vertices(self): return self.A.ordered_vertices() + self.B.ordered_vertices() @@ -1037,7 +1047,9 @@ def dim(self): def d_entities(self, d, get_class=True): if isinstance(d, tuple): - return [TensorProductPoint(e_a, e_b) for e_a in self.A.d_entities(d[0], get_class) for e_b in self.B.d_entities(d[1], get_class)] + if get_class: + return self.entities[d] + return [e.id for e in self.entities[d]] raise NotImplementedError("not sure this is right") return self.A.d_entities(d, get_class) + self.B.d_entities(d, get_class) diff --git a/fuse/dof.py b/fuse/dof.py index 4c6546f5..d6eab3a4 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -195,7 +195,7 @@ def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): else: comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] - if isinstance(self.pt, int): + if not isinstance(self.pt, tuple): return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), comps if not immersed: return Qpts, np.array([wt*np.matmul(self.pt, basis_change)for wt in Qwts]).astype(np.float64), comps diff --git a/fuse/triples.py b/fuse/triples.py index 4a68e544..4e7dc654 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -18,6 +18,20 @@ from functools import cache +def compute_form_degree(cell, spaces): + if str(spaces[1]) == "L2": + return cell.dimension + elif str(spaces[1]) == "H1": + return 0 + if cell.dimension < 2: + raise ValueError(f"Cells of dimension {cell.dimension} can only be 0 or 1 forms") + if cell.dimension == 2: + return 1 + elif str(spaces[1]) == "HDiv": + return 1 + elif str(spaces[1]) == "HCurl": + return 2 + class ElementTriple(): """ Class to represent the three core parts of the element @@ -157,7 +171,8 @@ def to_ufl(self): def to_fiat(self): # call this to ensure set up is complete self.to_ufl() - form_degree = 1 if self.spaces[0].set_shape else 0 + # form_degree = 1 if self.spaces[0].set_shape else 0 + form_degree = compute_form_degree(self.cell, self.spaces) degree = self.spaces[0].degree() # sanity check that the dofs span the space original_V, original_basis = self.compute_dense_matrix(self.ref_el, self.entity_ids, self.nodes, self.poly_set) @@ -544,8 +559,6 @@ def reverse_dof_perms(self, matrices): for e in ents: # old_e_id = e.id - min_ids[dim] e_id = self.cell.d_entities(e.dim(), get_class=False).index(e.id) - if isinstance(dim, tuple): - breakpoint() perms_copy = matrices[dim][e_id].copy() members = e.group.members() for m in members: diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index 1f9f1ca0..8cb17064 100644 --- a/test/test_2d_examples_docs.py +++ b/test/test_2d_examples_docs.py @@ -1,6 +1,7 @@ from fuse import * import sympy as sp import numpy as np +import pytest np.set_printoptions(legacy="1.25") @@ -32,7 +33,7 @@ def construct_dg1(): def construct_dg0_integral(cell=None): edge = Point(1, [Point(0), Point(0)], vertex_num=2) - xs = [DOF(L2Pairing(), VectorKernel(1))] + xs = [DOF(L2Pairing(), VectorKernel(0.5))] dg0 = ElementTriple(edge, (P0, CellL2, C0), DOFGenerator(xs, S1, S1)) return dg0 @@ -345,6 +346,15 @@ def test_rt_example(): rt.to_fiat() +@pytest.mark.parametrize(["triple", "expected"], [(construct_dg1_integral(), 1), + (construct_cg1(), 0), + (construct_rt(), 1), + (construct_dg1_tri(), 2)]) +def test_form_degree(triple, expected): + from fuse.triples import compute_form_degree + assert compute_form_degree(triple.cell, triple.spaces) == expected + + def construct_hermite(): tri = polygon(3) vert = tri.vertices()[0] diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 7fbb158b..fe2e24e6 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -2,7 +2,7 @@ import numpy as np from fuse import * from firedrake import * -from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg1_integral +from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg0_integral, construct_dg1_integral from test_convert_to_fiat import create_cg2 # from test_convert_to_fiat import create_cg1 @@ -148,6 +148,10 @@ def test_cg1_dg0(): # print(flatten_entities(non_sym.entity_dofs())) # print(flatten_entities(non_sym2.entity_dofs())) print(non_sym) + A = construct_dg1_integral() + B = construct_dg1_integral() + non_sym1 = tensor_product(A, B) + print(non_sym1) breakpoint() @@ -195,10 +199,16 @@ def test_hdiv(): mesh = ExtrudedMesh(m, 2) CG_1 = FiniteElement("CG", "interval", 1) DG_0 = FiniteElement("DG", "interval", 0) + cg1 = construct_cg1() + dg0 = construct_dg0_integral() + p1p0 = tensor_product(cg1, dg0).to_ufl() P1P0 = TensorProductElement(CG_1, DG_0) - RT_horiz = HDivElement(P1P0) + RT_horiz = HDivElement(p1p0) + # RT_horiz = HDivElement(P1P0) + p0p1 = tensor_product(dg0, cg1).to_ufl() P0P1 = TensorProductElement(DG_0, CG_1) - RT_vert = HDivElement(P0P1) + RT_vert = HDivElement(p0p1) + # RT_vert = HDivElement(P0P1) elt = RT_horiz + RT_vert # mesh = UnitSquareMesh(1, 1, quadrilateral=True) V = FunctionSpace(mesh, elt) @@ -206,7 +216,8 @@ def test_hdiv(): v = TestFunction(V) f = Function(V) x, y = SpatialCoordinate(mesh) - f_vec = as_vector(((1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2), (1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2))) + # f_vec = as_vector(((1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2), (1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2))) + f_vec = as_vector((2, 3)) f = project(f_vec, V) a = (inner(grad(u), grad(v)) + inner(u, v)) * dx L = inner(f, v) * dx From 9c107c0e198450458cb3682f3cd8da2ba276a7f6 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 4 May 2026 09:54:55 +0100 Subject: [PATCH 29/32] add higher order tests, turn matrix use off --- fuse/cells.py | 3 --- fuse/tensor_products.py | 22 +++++++++++------ fuse/triples.py | 3 ++- test/test_tensor_prod.py | 51 +++++++++++++++++++++++++++++++--------- 4 files changed, 57 insertions(+), 22 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index b7a9d1d4..62038eb4 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -986,9 +986,6 @@ def __init__(self, A, B): self.entities[d] = [self] else: self.entities[d] = [TensorProductPoint(e_a, e_b) for e_a in self.A.d_entities(d[0], True) for e_b in self.B.d_entities(d[1], True)] - - - def ordered_vertices(self): return self.A.ordered_vertices() + self.B.ordered_vertices() diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index 8891ce48..a37b7f69 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -4,15 +4,21 @@ from finat.ufl import TensorProductElement, FuseElement -def tensor_product(A, B): +def tensor_product(A, B, matrices=True): if not (isinstance(A, ElementTriple) and isinstance(B, ElementTriple)): raise ValueError("Both components of Tensor Product need to be a Fuse Triple.") - return TensorProductTriple(A, B) + return TensorProductTriple(A, B, matrices=matrices) + + +def symmetric_tensor_product(A, B): + if not (isinstance(A, ElementTriple) and isinstance(B, ElementTriple)): + raise ValueError("Both components of Tensor Product need to be a Fuse Triple.") + return TensorProductTriple(A, B, symmetric=True) class TensorProductTriple(ElementTriple): - def __init__(self, A, B, flat=False): + def __init__(self, A, B, flat=False, symmetric=True, matrices=True): self.A = A self.B = B self.spaces = [] @@ -21,12 +27,14 @@ def __init__(self, A, B, flat=False): self.DOFGenerator = [A.DOFGenerator, B.DOFGenerator] self.cell = TensorProductPoint(A.cell, B.cell) + self.symmetric = symmetric self.flat = flat - self.apply_matrices = False - self.setup_matrices() if self.flat: self.unflat_cell = self.cell self.cell = self.cell.flatten() + self.apply_matrices = matrices + if self.apply_matrices: + self.setup_matrices() def sub_elements(self): return [self.A, self.B] @@ -105,7 +113,7 @@ def to_ufl(self): return TensorProductElement(*ufl_sub_elements, cell=self.cell.to_ufl()) def flatten(self): - return TensorProductTriple(self.A, self.B, flat=True) + return TensorProductTriple(self.A, self.B, flat=True, symmetric=self.symmetric, matrices=self.apply_matrices) def unflatten(self): - return TensorProductTriple(self.A, self.B, flat=False) + return TensorProductTriple(self.A, self.B, flat=False, symmetric=self.symmetric, matrices=self.apply_matrices) diff --git a/fuse/triples.py b/fuse/triples.py index 4e7dc654..6f39aa4a 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -30,7 +30,8 @@ def compute_form_degree(cell, spaces): elif str(spaces[1]) == "HDiv": return 1 elif str(spaces[1]) == "HCurl": - return 2 + return 2 + class ElementTriple(): """ diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index fe2e24e6..1d5fafe7 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -3,10 +3,26 @@ from fuse import * from firedrake import * from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg0_integral, construct_dg1_integral -from test_convert_to_fiat import create_cg2 +from test_convert_to_fiat import create_cg2, create_dg0 # from test_convert_to_fiat import create_cg1 +def create_cg3_interval(cell=None): + if cell is None: + cell = line() + deg = 3 + if cell.dim() > 1: + raise NotImplementedError("This method is for cg3 on edges, please use construct_cg3 for triangles") + vert_dg = create_dg0(cell.vertices()[0]) + xs = [immerse(cell, vert_dg, TrH1)] + interior = [DOF(DeltaPairing(), PointKernel((-1/np.sqrt(5), )))] + + Pk = PolynomialSpace(deg) + cg = ElementTriple(cell, (Pk, CellL2, C0), [DOFGenerator(xs, get_cyc_group(len(cell.vertices())), S1), + DOFGenerator(interior, S2, S1)]) + return cg + + def helmholtz_solve(mesh, V): u = TrialFunction(V) v = TestFunction(V) @@ -63,15 +79,18 @@ def test_ext_mesh(generator1, generator2, code1, code2, deg1, deg2): assert np.allclose(res1, res2) -def test_helmholtz(): +@pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(construct_cg1, "CG", 1, 1.8), + (create_cg2, "CG", 2, 3.8), + (create_cg3_interval, "CG", 3, 4.8)]) +def test_helmholtz(elem_gen, elem_code, deg, conv_rate): vals = range(3, 6) res = [] for r in vals: m = UnitIntervalMesh(2**r) mesh = ExtrudedMesh(m, 2**r) - A = construct_cg1() - B = construct_cg1() + A = elem_gen() + B = elem_gen() elem = tensor_product(A, B) U = FunctionSpace(mesh, elem.to_ufl()) @@ -80,7 +99,7 @@ def test_helmholtz(): res = np.array(res) conv = np.log2(res[:-1] / res[1:]) print("convergence order:", conv) - assert (np.array(conv) > 1.8).all() + assert (np.array(conv) > conv_rate).all() def test_on_quad_mesh(): @@ -99,19 +118,23 @@ def test_on_quad_mesh(): @pytest.mark.parametrize(["elem_gen", "elem_code", "deg", "conv_rate"], [(construct_cg1, "CG", 1, 1.8), - (create_cg2, "CG", 2, 3.8)]) + (create_cg2, "CG", 2, 3.8), + (create_cg3_interval, "CG", 3, 4.8)]) def test_quad_mesh_helmholtz(elem_gen, elem_code, deg, conv_rate): quadrilateral = True vals = range(3, 6) res_fuse = [] + res_fire = [] for r in vals: mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) A = elem_gen() B = elem_gen() - elem = tensor_product(A, B).flatten() + elem = tensor_product(A, B, matrices=False).flatten() U = FunctionSpace(mesh, elem.to_ufl()) res_fuse += [helmholtz_solve(mesh, U)] + U = FunctionSpace(mesh, elem_code, deg) + res_fire += [helmholtz_solve(mesh, U)] print("Fuse l2 error norms:", res_fuse) res = np.array(res_fuse) @@ -119,6 +142,12 @@ def test_quad_mesh_helmholtz(elem_gen, elem_code, deg, conv_rate): print("Fuse convergence order:", conv) assert (np.array(conv) > conv_rate).all() + print("FIAT l2 error norms:", res_fire) + res = np.array(res_fire) + conv = np.log2(res[:-1] / res[1:]) + print("Fiat convergence order:", conv) + assert (np.array(conv) > conv_rate).all() + @pytest.mark.parametrize(["A", "B", "res"], [(Point(0), line(), False), (line(), line(), True), @@ -197,16 +226,16 @@ def test_hdiv(): np.set_printoptions(linewidth=90, precision=4, suppress=True) m = UnitIntervalMesh(2) mesh = ExtrudedMesh(m, 2) - CG_1 = FiniteElement("CG", "interval", 1) - DG_0 = FiniteElement("DG", "interval", 0) + # CG_1 = FiniteElement("CG", "interval", 1) + # DG_0 = FiniteElement("DG", "interval", 0) cg1 = construct_cg1() dg0 = construct_dg0_integral() p1p0 = tensor_product(cg1, dg0).to_ufl() - P1P0 = TensorProductElement(CG_1, DG_0) + # P1P0 = TensorProductElement(CG_1, DG_0) RT_horiz = HDivElement(p1p0) # RT_horiz = HDivElement(P1P0) p0p1 = tensor_product(dg0, cg1).to_ufl() - P0P1 = TensorProductElement(DG_0, CG_1) + # P0P1 = TensorProductElement(DG_0, CG_1) RT_vert = HDivElement(p0p1) # RT_vert = HDivElement(P0P1) elt = RT_horiz + RT_vert From 55ac2f9b07c43d90d2cd54f3d7730ae9089a40a6 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 4 May 2026 16:46:39 +0100 Subject: [PATCH 30/32] adapt tensor product creation to follow fuse directions, include in matrices --- fuse/tensor_products.py | 138 ++++++++++++++++++++++++++++++++-- test/test_2d_examples_docs.py | 11 ++- test/test_tensor_prod.py | 45 +++++++++-- 3 files changed, 177 insertions(+), 17 deletions(-) diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index a37b7f69..6113d718 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -1,7 +1,8 @@ -from fuse.triples import ElementTriple +from fuse.triples import ElementTriple, compute_form_degree +from fuse.traces import TrHCurl, TrHDiv from fuse.cells import TensorProductPoint import numpy as np -from finat.ufl import TensorProductElement, FuseElement +from finat.ufl import TensorProductElement, FuseElement, HDivElement, HCurlElement def tensor_product(A, B, matrices=True): @@ -16,6 +17,12 @@ def symmetric_tensor_product(A, B): return TensorProductTriple(A, B, symmetric=True) +def hcurl_transform(tensor_element): + gem_transformer, mat_transformer = select_fuse_hcurl_transformer(tensor_element) + tensor_element.add_mat_transformer(mat_transformer, TrHCurl) + return gem_transformer + + class TensorProductTriple(ElementTriple): def __init__(self, A, B, flat=False, symmetric=True, matrices=True): @@ -32,10 +39,12 @@ def __init__(self, A, B, flat=False, symmetric=True, matrices=True): if self.flat: self.unflat_cell = self.cell self.cell = self.cell.flatten() + self.mat_transformer = None self.apply_matrices = matrices if self.apply_matrices: self.setup_matrices() + @property def sub_elements(self): return [self.A, self.B] @@ -61,7 +70,7 @@ def setup_matrices(self): comp_os = self.cell.component_orientations() for e, (a, b) in enumerate(ents): ent_dofs = self.entity_dofs[dim][(a, b)] - if len(ent_dofs) > 1: + if len(ent_dofs) >= 1: sub_mat = oriented_mats_by_entity[dim][e] a_mat = self.A.matrices[dim[0]][a] a_ent_ids = self.A.entity_ids[dim[0]][a] @@ -72,7 +81,11 @@ def setup_matrices(self): for o in os: a_sub_mat = a_mat[o[0]][np.ix_(a_ent_ids, a_ent_ids)] b_sub_mat = b_mat[o[1]][np.ix_(b_ent_ids, b_ent_ids)] - combined_sub_mat = np.kron(a_sub_mat, b_sub_mat) + if self.mat_transformer is not None: + o_classes = (self.A.cell.group.get_member_by_val(o[0]), self.B.cell.group.get_member_by_val(o[1])) + combined_sub_mat = self.mat_transformer(a_sub_mat, b_sub_mat, o_classes) + else: + combined_sub_mat = np.kron(a_sub_mat, b_sub_mat) new_o = comp_os[dim][o] sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)] = np.matmul(sub_mat[new_o][np.ix_(ent_dofs, ent_dofs)], combined_sub_mat) # from collections import defaultdict @@ -109,11 +122,124 @@ def generate(self): def to_ufl(self): if self.flat: return FuseElement(self, self.cell.to_ufl()) - ufl_sub_elements = [e.to_ufl() for e in self.sub_elements()] - return TensorProductElement(*ufl_sub_elements, cell=self.cell.to_ufl()) + ufl_sub_elements = [e.to_ufl() for e in self.sub_elements] + return TensorProductElement(*ufl_sub_elements, cell=self.cell.to_ufl(), triple=self) def flatten(self): return TensorProductTriple(self.A, self.B, flat=True, symmetric=self.symmetric, matrices=self.apply_matrices) def unflatten(self): return TensorProductTriple(self.A, self.B, flat=False, symmetric=self.symmetric, matrices=self.apply_matrices) + + +def compute_matrix_transform(trace, cell, o): + bvs = np.array(cell.basis_vectors()) + new_bvs = np.array(cell.orient(~o).basis_vectors()) + basis_change = np.matmul(new_bvs, np.linalg.inv(bvs)) + # if len(ent_dofs_ids) == basis_change.shape[0]: + # sub_mat = basis_change + # elif len(dof_gen_class[dim].g2.members()) == 2 and len(ent_dofs_ids) == 1: + # # equivalently g1 trivial + # sub_mat = trace.manipulate_basis(basis_change) + # else: + # case where value change is a restriction of the full transformation of the basis + value_change = trace(cell).manipulate_basis(basis_change) + # sub_mat = np.kron((~o).matrix_form(), value_change) + return value_change + + +class HDiv(TensorProductTriple): + + def __init__(self, tensor_element): + self.gem_transformer, self.mat_transformer = self.select_fuse_hdiv_transformer(tensor_element) + self.trace = TrHDiv + super(HDiv, self).__init__(tensor_element.A, tensor_element.B, tensor_element.flat, tensor_element.symmetric, tensor_element.matrices) + self.spaces[1] = TrHDiv + + def to_ufl(self): + return HDivElement(super(HDiv, self).to_ufl(), self.gem_transformer) + + def select_fuse_hdiv_transformer(self, element): + # Assume: something x interval + import gem + assert len(element.sub_elements) == 2 + assert element.sub_elements[1].cell.get_shape() == 1 + # Globally consistent edge orientations of the reference + # quadrilateral: rightward horizontally, upward vertically. + # Their rotation by 90 degrees anticlockwise is interpreted as the + # positive direction for normal vectors. + ks = tuple(compute_form_degree(fe.cell, fe.spaces) for fe in element.sub_elements) + transform = lambda cell, o: compute_matrix_transform(element.trace, cell, o) + if ks == (0, 1): + # Make the scalar value the right hand rule normal on the + # y-aligned edges. + cell = element.sub_elements[1].cell + bv = cell.basis_vectors()[0] + return lambda v: [gem.Product(gem.Literal(bv[0]), v), gem.Zero()], lambda m_a, m_b, o: np.kron(transform(cell, o[1]) @ m_a, m_b) + elif ks == (1, 0): + # Make the scalar value the upward-pointing normal on the + # x-aligned edges. + cell = element.sub_elements[0].cell + bv = cell.basis_vectors()[0][0] + return lambda v: [gem.Zero(), gem.Product(gem.Literal(bv), v)], lambda m_a, m_b, o: np.kron(m_a, transform(cell, o[0]) @ m_b) + # elif ks == (2, 0): + # # Same for 3D, so z-plane. + # return lambda v: [gem.Zero(), gem.Zero(), v] + # elif ks == (1, 1): + # if element.mapping == "contravariant piola": + # # Pad the 2-vector normal on the "base" cell into a + # # 3-vector, maintaining direction. + # return lambda v: [gem.Indexed(v, (0,)), + # gem.Indexed(v, (1,)), + # gem.Zero()] + # elif element.mapping == "covariant piola": + # # Rotate the 2-vector tangential component on the "base" + # # cell 90 degrees anticlockwise into a 3-vector and pad. + # return lambda v: [gem.Indexed(v, (1,)), + # gem.Product(gem.Literal(-1), gem.Indexed(v, (0,))), + # gem.Zero()] + # else: + # assert False, "Unexpected original mapping!" + else: + raise NotImplementedError("Unexpected original mapping!") + assert False, "Unexpected form degree combination!" + + def select_fuse_hcurl_transformer(element): + import gem + # Assume: something x interval + assert len(element.sub_elements) == 2 + assert element.sub_elements[1].cell.get_shape() == 1 + + # Globally consistent edge orientations of the reference + # quadrilateral: rightward horizontally, upward vertically. + # Tangential vectors interpret these as the positive direction. + dim = element.cell.get_spatial_dimension() + ks = tuple(compute_form_degree(fe.cell, fe.spaces) for fe in element.sub_elements) + if all(str(fe.spaces[1]) == "H1" or str(fe.spaces[1]) == "L2" for fe in element.sub_elements): # affine mapping + if ks == (1, 0): + # Can only be 2D. Make the scalar value the + # tangential following the cell edge direction on the x-aligned edges. + bv = element.sub_elements[0].cell.basis_vectors()[0][0] + return lambda v: [gem.Product(gem.Literal(bv), v), gem.Zero()] + elif ks == (0, 1): + # Can be any spatial dimension. Make the scalar value the + # tangential following the cell edge direction . + bv = element.sub_elements[1].cell.basis_vectors()[0][0] + return lambda v: [gem.Zero()] * (dim - 1) + [gem.Product(gem.Literal(bv), v)] + else: + assert False + # elif any(str(fe.spaces[1]) == "HCurl" for fe in element.sub_elements): # Covariant Piola mapping + # # Second factor must be continuous interval. Just padding. + # return lambda v: [gem.Indexed(v, (0,)), + # gem.Indexed(v, (1,)), + # gem.Zero()] + # elif any(str(fe.spaces[1]) == "HDiv" for fe in element.sub_elements): # Contravariant Piola mapping + # # Second factor must be continuous interval. Rotate the + # # 2-vector tangential component on the "base" cell 90 degrees + # # clockwise into a 3-vector and pad. + # return lambda v: [gem.Product(gem.Literal(-1), gem.Indexed(v, (1,))), + # gem.Indexed(v, (0,)), + # gem.Zero()] + else: + raise NotImplementedError("Unexpected original mapping!") + assert False, "Unexpected original mapping!" diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index deec1a48..9754bb42 100644 --- a/test/test_2d_examples_docs.py +++ b/test/test_2d_examples_docs.py @@ -31,8 +31,9 @@ def construct_dg1(): return dg1 -def construct_dg0_integral(cell=None): - edge = Point(1, [Point(0), Point(0)], vertex_num=2) +def construct_dg0_integral(edge=None): + if not edge: + edge = Point(1, [Point(0), Point(0)], vertex_num=2) xs = [DOF(L2Pairing(), VectorKernel(0.5))] dg0 = ElementTriple(edge, (P0, CellL2, C0), DOFGenerator(xs, S1, S1)) return dg0 @@ -105,9 +106,11 @@ def test_dg_examples(): assert any(np.isclose(val, dof.eval(test_func)) for val in dof_vals) -def construct_cg1(): +def construct_cg1(edge=None): + # [test_cg1 0] - edge = Point(1, [Point(0), Point(0)], vertex_num=2) + if not edge: + edge = Point(1, [Point(0), Point(0)], vertex_num=2) vert = edge.vertices()[0] xs = [DOF(DeltaPairing(), PointKernel(()))] diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 1d5fafe7..1c15ea3f 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -223,22 +223,26 @@ def test_trace_galerkin_projection(): def test_hdiv(): + from fuse.tensor_products import HDiv np.set_printoptions(linewidth=90, precision=4, suppress=True) m = UnitIntervalMesh(2) mesh = ExtrudedMesh(m, 2) - # CG_1 = FiniteElement("CG", "interval", 1) - # DG_0 = FiniteElement("DG", "interval", 0) + CG_1 = FiniteElement("CG", "interval", 1) + DG_0 = FiniteElement("DG", "interval", 0) cg1 = construct_cg1() dg0 = construct_dg0_integral() - p1p0 = tensor_product(cg1, dg0).to_ufl() + p1p0 = HDiv(tensor_product(cg1, dg0)) # P1P0 = TensorProductElement(CG_1, DG_0) - RT_horiz = HDivElement(p1p0) + # RT_horiz = HDivElement(p1p0.to_ufl(), transform=hdiv_transform(p1p0)) + RT_horiz = p1p0.to_ufl() # RT_horiz = HDivElement(P1P0) - p0p1 = tensor_product(dg0, cg1).to_ufl() + p0p1 = HDiv(tensor_product(dg0, cg1)) # P0P1 = TensorProductElement(DG_0, CG_1) - RT_vert = HDivElement(p0p1) + RT_vert = p0p1.to_ufl() # RT_vert = HDivElement(P0P1) - elt = RT_horiz + RT_vert + elt = RT_horiz + # + RT_vert + # + RT_vert # mesh = UnitSquareMesh(1, 1, quadrilateral=True) V = FunctionSpace(mesh, elt) u = TrialFunction(V) @@ -260,3 +264,30 @@ def test_hdiv(): # print(ent) # for comp in arr: # print(comp[0], comp[1]) + + +def test_transforms(): + edge = Point(1, [Point(0), Point(0)], vertex_num=2) + rev_edge = edge.orient(edge.group.members()[1]) + from fuse.tensor_products import hcurl_transform, hdiv_transform + cg1 = construct_cg1() + rev_cg1 = construct_cg1(rev_edge) + dg0 = construct_dg0_integral() + rev_dg0 = construct_dg0_integral(rev_edge) + import gem + v = gem.Literal(5) + # print("HCurl") + # print(hcurl_transform(tensor_product(dg0, cg1))(v)) + # print(hcurl_transform(tensor_product(rev_dg0, cg1))(v)) + # print(hcurl_transform(tensor_product(cg1, dg0))(v)) + # print(hcurl_transform(tensor_product(cg1, rev_dg0))(v)) + # print(hcurl_transform(tensor_product(dg0, rev_cg1))(v)) + # print(hcurl_transform(tensor_product(rev_cg1, dg0))(v)) + print("HDiv") + print(hdiv_transform(tensor_product(dg0, cg1))(v)) + print(hdiv_transform(tensor_product(rev_dg0, cg1))(v)) + print(hdiv_transform(tensor_product(cg1, dg0))(v)) + print(hdiv_transform(tensor_product(cg1, rev_dg0))(v)) + print(hdiv_transform(tensor_product(dg0, rev_cg1))(v)) + print(hdiv_transform(tensor_product(rev_cg1, dg0))(v)) + breakpoint() From 5748785a952660b193609a628a42f4b3f29f5045 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 4 May 2026 16:50:45 +0100 Subject: [PATCH 31/32] lint --- fuse/tensor_products.py | 18 ++++++++++++------ test/test_perms.py | 2 +- test/test_tensor_prod.py | 38 +++++++++++++++++++------------------- 3 files changed, 32 insertions(+), 26 deletions(-) diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index 6113d718..25407856 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -17,12 +17,6 @@ def symmetric_tensor_product(A, B): return TensorProductTriple(A, B, symmetric=True) -def hcurl_transform(tensor_element): - gem_transformer, mat_transformer = select_fuse_hcurl_transformer(tensor_element) - tensor_element.add_mat_transformer(mat_transformer, TrHCurl) - return gem_transformer - - class TensorProductTriple(ElementTriple): def __init__(self, A, B, flat=False, symmetric=True, matrices=True): @@ -204,6 +198,18 @@ def select_fuse_hdiv_transformer(self, element): raise NotImplementedError("Unexpected original mapping!") assert False, "Unexpected form degree combination!" + +class HCurl(TensorProductTriple): + + def __init__(self, tensor_element): + self.gem_transformer, self.mat_transformer = self.select_fuse_hcurl_transformer(tensor_element) + self.trace = TrHCurl + super(HDiv, self).__init__(tensor_element.A, tensor_element.B, tensor_element.flat, tensor_element.symmetric, tensor_element.matrices) + self.spaces[1] = TrHCurl + + def to_ufl(self): + return HCurlElement(super(HDiv, self).to_ufl(), self.gem_transformer) + def select_fuse_hcurl_transformer(element): import gem # Assume: something x interval diff --git a/test/test_perms.py b/test/test_perms.py index ec403774..25c71231 100644 --- a/test/test_perms.py +++ b/test/test_perms.py @@ -1,5 +1,5 @@ from fuse import * -from test_convert_to_fiat import create_cg1, create_dg1, create_cg2, construct_nd +from test_convert_to_fiat import create_cg1, create_dg1, create_cg2 from test_2d_examples_docs import construct_cg3 import pytest import numpy as np diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 1c15ea3f..713c9621 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -227,8 +227,8 @@ def test_hdiv(): np.set_printoptions(linewidth=90, precision=4, suppress=True) m = UnitIntervalMesh(2) mesh = ExtrudedMesh(m, 2) - CG_1 = FiniteElement("CG", "interval", 1) - DG_0 = FiniteElement("DG", "interval", 0) + # CG_1 = FiniteElement("CG", "interval", 1) + # DG_0 = FiniteElement("DG", "interval", 0) cg1 = construct_cg1() dg0 = construct_dg0_integral() p1p0 = HDiv(tensor_product(cg1, dg0)) @@ -236,11 +236,11 @@ def test_hdiv(): # RT_horiz = HDivElement(p1p0.to_ufl(), transform=hdiv_transform(p1p0)) RT_horiz = p1p0.to_ufl() # RT_horiz = HDivElement(P1P0) - p0p1 = HDiv(tensor_product(dg0, cg1)) + # p0p1 = HDiv(tensor_product(dg0, cg1)) # P0P1 = TensorProductElement(DG_0, CG_1) - RT_vert = p0p1.to_ufl() + # RT_vert = p0p1.to_ufl() # RT_vert = HDivElement(P0P1) - elt = RT_horiz + elt = RT_horiz # + RT_vert # + RT_vert # mesh = UnitSquareMesh(1, 1, quadrilateral=True) @@ -269,25 +269,25 @@ def test_hdiv(): def test_transforms(): edge = Point(1, [Point(0), Point(0)], vertex_num=2) rev_edge = edge.orient(edge.group.members()[1]) - from fuse.tensor_products import hcurl_transform, hdiv_transform + from fuse.tensor_products import HDiv, HCurl cg1 = construct_cg1() rev_cg1 = construct_cg1(rev_edge) dg0 = construct_dg0_integral() rev_dg0 = construct_dg0_integral(rev_edge) import gem v = gem.Literal(5) - # print("HCurl") - # print(hcurl_transform(tensor_product(dg0, cg1))(v)) - # print(hcurl_transform(tensor_product(rev_dg0, cg1))(v)) - # print(hcurl_transform(tensor_product(cg1, dg0))(v)) - # print(hcurl_transform(tensor_product(cg1, rev_dg0))(v)) - # print(hcurl_transform(tensor_product(dg0, rev_cg1))(v)) - # print(hcurl_transform(tensor_product(rev_cg1, dg0))(v)) + print("HCurl") + print(HCurl(tensor_product(dg0, cg1))(v)) + print(HCurl(tensor_product(rev_dg0, cg1))(v)) + print(HCurl(tensor_product(cg1, dg0))(v)) + print(HCurl(tensor_product(cg1, rev_dg0))(v)) + print(HCurl(tensor_product(dg0, rev_cg1))(v)) + print(HCurl(tensor_product(rev_cg1, dg0))(v)) print("HDiv") - print(hdiv_transform(tensor_product(dg0, cg1))(v)) - print(hdiv_transform(tensor_product(rev_dg0, cg1))(v)) - print(hdiv_transform(tensor_product(cg1, dg0))(v)) - print(hdiv_transform(tensor_product(cg1, rev_dg0))(v)) - print(hdiv_transform(tensor_product(dg0, rev_cg1))(v)) - print(hdiv_transform(tensor_product(rev_cg1, dg0))(v)) + print(HDiv(tensor_product(dg0, cg1))(v)) + print(HDiv(tensor_product(rev_dg0, cg1))(v)) + print(HDiv(tensor_product(cg1, dg0))(v)) + print(HDiv(tensor_product(cg1, rev_dg0))(v)) + print(HDiv(tensor_product(dg0, rev_cg1))(v)) + print(HDiv(tensor_product(rev_cg1, dg0))(v)) breakpoint() From 7d03fb51d88c4dc46a2a1e21f6541c180ba7f3b7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 5 May 2026 12:49:27 +0100 Subject: [PATCH 32/32] fix merge issues --- fuse/dof.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fuse/dof.py b/fuse/dof.py index 83b90ffd..2c6cad0a 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -195,8 +195,8 @@ def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): else: comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] - if isinstance(self.pt, int): - return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] + if not isinstance(self.pt, tuple): + return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), comps if not immersed: return Qpts, np.array([wt*np.matmul(self.pt, basis_change)for wt in Qwts]).astype(np.float64), comps return Qpts, np.array([wt*immersed(np.matmul(self.pt, basis_change))for wt in Qwts]).astype(np.float64), comps