From 146ac89123514ac27d88dd6c6bab7f0970b41e85 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 17 Jan 2025 16:29:12 +0000 Subject: [PATCH 01/10] Changes to fix convergence of firedrake CG... now to fix mine --- fuse/cells.py | 8 +++++++- fuse/triples.py | 5 +++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index c188de0..b752c58 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -383,13 +383,17 @@ def get_topology(self): self.topology = {} self.topology_unrelabelled = {} + self.topology_unrelabelled = {} for i in range(len(structure)): dimension = structure[i] self.topology[i] = {} self.topology_unrelabelled[i] = {} + self.topology_unrelabelled[i] = {} for node in dimension: self.topology[i][node - min_ids[i]] = tuple([relabelled_verts[vert] for vert in self.get_node(node).ordered_vertices()]) self.topology_unrelabelled[i][node - min_ids[i]] = tuple([vert - min_ids[0] for vert in self.get_node(node).ordered_vertices()]) + return self.topology_unrelabelled + self.topology_unrelabelled[i][node - min_ids[i]] = tuple([vert - min_ids[0] for vert in self.get_node(node).ordered_vertices()]) return self.topology_unrelabelled def get_starter_ids(self): @@ -794,6 +798,8 @@ def __init__(self, cell, name=None): name = "FuseCell" self.name = name + # vertices = sorted(cell.ordered_vertices()) + # verts = [cell.get_node(v, return_coords=True) for v in cell.ordered_vertices()] verts = cell.vertices(return_coords=True) topology = cell.get_topology() @@ -936,7 +942,7 @@ def to_fiat(self): return self.cell_complex.to_fiat(name=self.cellname()) def __repr__(self): - return super(CellComplexToUFL, self).__repr__() + return super(CellComplexToUFL, self).__repr__() def reconstruct(self, **kwargs): """Reconstruct this cell, overwriting properties by those in kwargs.""" diff --git a/fuse/triples.py b/fuse/triples.py index f977b2a..b2b9cc0 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -114,12 +114,13 @@ def to_fiat(self): entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])] # if sort_entities: # # sort the entities by support vertex ids - # support = [sorted(top[dim][entity]) for dim, entity in entities] - # entities = [entity for verts, entity in sorted(zip(support, entities))] + support = [sorted(top[dim][entity]) for dim, entity in entities] + entities = [entity for verts, entity in sorted(zip(support, entities))] counter = 0 for entity in entities: dim = entity[0] for i in range(len(dofs)): + # breakpoint() if entity[1] == dofs[i].trace_entity.id - min_ids[dim]: entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].append(counter) nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) From 9e8da69644f5ccf82387bf22d94e52909c350c08 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 20 Jan 2025 15:36:06 +0000 Subject: [PATCH 02/10] Debugging issues with lagrange --- fuse/cells.py | 2 +- fuse/triples.py | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index b752c58..33133d0 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -798,8 +798,8 @@ def __init__(self, cell, name=None): name = "FuseCell" self.name = name - # vertices = sorted(cell.ordered_vertices()) + # verts = [cell.get_node(v, return_coords=True) for v in cell.ordered_vertices()] verts = cell.vertices(return_coords=True) topology = cell.get_topology() diff --git a/fuse/triples.py b/fuse/triples.py index b2b9cc0..f977b2a 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -114,13 +114,12 @@ def to_fiat(self): entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])] # if sort_entities: # # sort the entities by support vertex ids - support = [sorted(top[dim][entity]) for dim, entity in entities] - entities = [entity for verts, entity in sorted(zip(support, entities))] + # support = [sorted(top[dim][entity]) for dim, entity in entities] + # entities = [entity for verts, entity in sorted(zip(support, entities))] counter = 0 for entity in entities: dim = entity[0] for i in range(len(dofs)): - # breakpoint() if entity[1] == dofs[i].trace_entity.id - min_ids[dim]: entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].append(counter) nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) From c6289d907206aea4101c2a1678951a936f8895c1 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 27 Jan 2025 10:44:54 +0000 Subject: [PATCH 03/10] Demo that flipping cells breaks CG3, worked out mocking, some other test stuff --- Makefile | 7 +++++++ fuse/cells.py | 1 + 2 files changed, 8 insertions(+) diff --git a/Makefile b/Makefile index b09392e..4360598 100644 --- a/Makefile +++ b/Makefile @@ -36,4 +36,11 @@ test_cells: @firedrake-clean @python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect1] +test_cells: + @echo " Running all cell comparison tests" + @firedrake-clean + @python -m pytest -rPx test/test_cells.py::test_ref_els[expect0] + @firedrake-clean + @python -m pytest -rPx test/test_cells.py::test_ref_els[expect1] + prepush: lint tests doc \ No newline at end of file diff --git a/fuse/cells.py b/fuse/cells.py index 33133d0..b487606 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -980,3 +980,4 @@ def constructCellComplex(name): return TensorProductPoint(*components).to_ufl(name) else: raise TypeError("Cell complex construction undefined for {}".format(str(name))) + From 53a52960c8510107eee69eb438a6e55094d536b4 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 27 Jan 2025 17:59:56 +0000 Subject: [PATCH 04/10] first thoughts about derivative/tangential fiat nodes --- fuse/dof.py | 28 +++++++++++++++++--- test/test_convert_to_fiat.py | 50 ++++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+), 4 deletions(-) diff --git a/fuse/dof.py b/fuse/dof.py index 4de0e66..9065282 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -37,6 +37,15 @@ def __call__(self, kernel, v, cell): def convert_to_fiat(self, ref_el, dof, interpolant_deg): pt = dof.eval(MyTestFunction(lambda *x: x)) return PointEvaluation(ref_el, pt) + + def get_pts(self, ref_el, total_degree): + entity = ref_el.construct_subelement(self.entity.dim()) + Q_ref = create_quadrature(entity, total_deg) + + ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] + Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref) + qpts, qwts = Q.get_points(), Q.get_weights() + return [(1,)], [(1,)] def add_entity(self, entity): res = DeltaPairing() @@ -90,6 +99,15 @@ def add_entity(self, entity): res.entity = entity return res + def get_pts(self, ref_el, total_degree): + entity = ref_el.construct_subelement(self.entity.dim()) + Q_ref = create_quadrature(entity, total_deg) + + ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] + Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref) + qpts, qwts = Q.get_points(), Q.get_weights() + return qpts, qwts + def convert_to_fiat(self, ref_el, dof, interpolant_degree): total_deg = interpolant_degree + dof.kernel.degree() ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] @@ -98,11 +116,7 @@ def convert_to_fiat(self, ref_el, dof, interpolant_degree): Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref) Jdet = Q.jacobian_determinant() qpts, _ = Q.get_points(), Q.get_weights() - print(qpts) - print(dof.tabulate(qpts)) f_at_qpts = dof.tabulate(qpts).T / Jdet - print(len(Q.pts)) - print(f_at_qpts.shape) functional = FrobeniusIntegralMoment(ref_el, Q, f_at_qpts) return functional @@ -295,7 +309,13 @@ def eval(self, fn, pullback=True): def tabulate(self, Qpts): immersion = self.target_space.tabulate(Qpts, self.trace_entity, self.g) res = self.kernel.tabulate(Qpts) + # [self.attachment(*tuple(r)) for r in res] return immersion*res + + def convert_to_fiat(self, ref_el): + pts, wts = self.pairing.get_pts() + self.target_space.convert_to_fiat(tabulated, wts) + def __call__(self, g): permuted = self.cell.permute_entities(g, self.trace_entity.dim()) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 4c00cf7..ee05805 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -154,6 +154,28 @@ def create_cg2_tri(cell): DOFGenerator(edge_xs, C3, S1)]) return cg +def create_hermite(tri): + vert = tri.vertices()[0] + + xs = [DOF(DeltaPairing(), PointKernel(()))] + dg0 = ElementTriple(vert, (P0, CellL2, C0), DOFGenerator(xs, S1, S1)) + + v_xs = [immerse(tri, dg0, TrH1)] + v_dofs = DOFGenerator(v_xs, S3/S2, S1) + + v_derv_xs = [immerse(tri, dg0, TrGrad)] + v_derv_dofs = DOFGenerator(v_derv_xs, S3/S2, S1) + + v_derv2_xs = [immerse(tri, dg0, TrHess)] + v_derv2_dofs = DOFGenerator(v_derv2_xs, S3/S2, S1) + + i_xs = [DOF(DeltaPairing(), PointKernel((0, 0)))] + i_dofs = DOFGenerator(i_xs, S1, S1) + + her = ElementTriple(tri, (P3, CellH2, C0), + [v_dofs, v_derv_dofs, v_derv2_dofs, i_dofs]) + return her + def test_create_fiat_nd(): cell = polygon(3) @@ -557,6 +579,34 @@ def test_project_3d(elem_gen, elem_code, deg): assert np.allclose(out.dat.data, f.dat.data, rtol=1e-5) +def test_create_hermite(): + deg = 3 + cell = polygon(3) + elem = create_hermite(cell) + ref_el = cell.to_fiat() + sd = ref_el.get_spatial_dimension() + + from FIAT.hermite import CubicHermite + fiat_elem = CubicHermite(ref_el, deg) + + my_elem = elem.to_fiat() + + Q = create_quadrature(ref_el, 2*(deg+1)) + Qpts, _ = Q.get_points(), Q.get_weights() + + fiat_vals = fiat_elem.tabulate(0, Qpts) + # my_vals = my_elem.tabulate(0, Qpts) + + fiat_vals = flatten(fiat_vals[(0,) * sd]) + # my_vals = flatten(my_vals[(0,) * sd]) + + # (x, res, _, _) = np.linalg.lstsq(fiat_vals.T, my_vals.T) + # x1 = np.linalg.inv(x) + # assert np.allclose(np.linalg.norm(my_vals.T - fiat_vals.T @ x), 0) + # assert np.allclose(np.linalg.norm(fiat_vals.T - my_vals.T @ x1), 0) + # assert np.allclose(res, 0) + + def test_investigate_dpc(): mesh = UnitSquareMesh(2, 2, quadrilateral=True) From 04271297c98562528cb678bfd0cc6d66b50f7416 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 28 Jan 2025 15:49:32 +0000 Subject: [PATCH 05/10] Refactor to allow derivative dof conversion --- fuse/dof.py | 69 ++++++++++++++++++++++++++---------- fuse/traces.py | 21 +++++++++++ fuse/triples.py | 16 +++++++-- test/test_convert_to_fiat.py | 2 +- 4 files changed, 86 insertions(+), 22 deletions(-) diff --git a/fuse/dof.py b/fuse/dof.py index 9065282..65e408f 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -1,6 +1,6 @@ from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule -from FIAT.functional import PointEvaluation, FrobeniusIntegralMoment +from FIAT.functional import PointEvaluation, FrobeniusIntegralMoment, Functional from fuse.utils import sympy_to_numpy import numpy as np import sympy as sp @@ -40,12 +40,7 @@ def convert_to_fiat(self, ref_el, dof, interpolant_deg): def get_pts(self, ref_el, total_degree): entity = ref_el.construct_subelement(self.entity.dim()) - Q_ref = create_quadrature(entity, total_deg) - - ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] - Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref) - qpts, qwts = Q.get_points(), Q.get_weights() - return [(1,)], [(1,)] + return [(0,) * entity.get_spatial_dimension()], [1], 1 def add_entity(self, entity): res = DeltaPairing() @@ -101,12 +96,15 @@ def add_entity(self, entity): def get_pts(self, ref_el, total_degree): entity = ref_el.construct_subelement(self.entity.dim()) - Q_ref = create_quadrature(entity, total_deg) + Q_ref = create_quadrature(entity, total_degree) ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref) - qpts, qwts = Q.get_points(), Q.get_weights() - return qpts, qwts + Jdet = Q.jacobian_determinant() + # TODO work out how to get J from attachment + + qpts, qwts = Q_ref.get_points(), Q_ref.get_weights() + return qpts, qwts, Jdet def convert_to_fiat(self, ref_el, dof, interpolant_degree): total_deg = interpolant_degree + dof.kernel.degree() @@ -168,8 +166,11 @@ def permute(self, g): def __call__(self, *args): return self.pt - def tabulate(self, Qpts): - return np.array([self.pt for _ in Qpts]).astype(np.float64) + def tabulate(self, Qpts, attachment=None): + if attachment is None: + return np.array([self.pt for _ in Qpts]).astype(np.float64) + else: + return np.array([attachment(*self.pt) for _ in Qpts]).astype(np.float64) def _to_dict(self): o_dict = {"pt": self.pt} @@ -209,8 +210,10 @@ def __call__(self, *args): return [res] return res - def tabulate(self, Qpts): - return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def tabulate(self, Qpts, attachment=None): + if attachment is None: + return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + return np.array([self(*attachment(*pt)) for pt in Qpts]).astype(np.float64) def _to_dict(self): o_dict = {"fn": self.fn} @@ -270,7 +273,23 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non def convert_to_fiat(self, ref_el, interpolant_degree): return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree) - raise NotImplementedError("Fiat conversion only implemented for Point eval") + + def convert_to_fiat_new(self, ref_el, interpolant_degree): + total_degree = self.kernel.degree() + interpolant_degree + pts, wts, jdet = self.pairing.get_pts(ref_el, total_degree) + f_pts = self.kernel.tabulate(pts).T / jdet + # TODO need to work out how i can discover the shape in a better way + if isinstance(self.pairing, DeltaPairing): + shp = tuple() + pt_dict = {tuple(p) : [(w, tuple())] for (p, w) in zip(f_pts.T, wts)} + else: + shp = tuple(f_pts.shape[:-1]) + weights = np.transpose(np.multiply(f_pts, wts), (-1,) + tuple(range(len(shp)))) + alphas = list(np.ndindex(shp)) + pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(pts, weights)} + + return Functional(ref_el, shp, pt_dict, {}, self.__repr__()) + def __repr__(self, fn="v"): return str(self.pairing).format(fn=fn, kernel=self.kernel) @@ -309,13 +328,25 @@ def eval(self, fn, pullback=True): def tabulate(self, Qpts): immersion = self.target_space.tabulate(Qpts, self.trace_entity, self.g) res = self.kernel.tabulate(Qpts) - # [self.attachment(*tuple(r)) for r in res] return immersion*res - def convert_to_fiat(self, ref_el): - pts, wts = self.pairing.get_pts() - self.target_space.convert_to_fiat(tabulated, wts) + def convert_to_fiat_new(self, ref_el, interpolant_degree): + total_degree = self.kernel.degree() + interpolant_degree + pts, wts, jdet = self.pairing.get_pts(ref_el, total_degree) + f_pts = self.kernel.tabulate(pts, self.attachment) + attached_pts = [self.attachment(*p) for p in pts] + immersion = self.target_space.tabulate(f_pts, self.trace_entity, self.g) + f_pts = (f_pts * immersion).T / jdet + pt_dict, deriv_dict = self.target_space.convert_to_fiat(attached_pts, f_pts, wts) + + # breakpoint() + # TODO need to work out how i can discover the shape in a better way + if isinstance(self.pairing, DeltaPairing): + shp = tuple() + else: + shp = tuple(f_pts.shape[:-1]) + return Functional(ref_el, shp, pt_dict, deriv_dict, self.__repr__()) def __call__(self, g): permuted = self.cell.permute_entities(g, self.trace_entity.dim()) diff --git a/fuse/traces.py b/fuse/traces.py index 46107ee..3d9f07b 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -54,6 +54,10 @@ def plot(self, ax, coord, trace_entity, g, **kwargs): def tabulate(self, Qpts, trace_entity, g): return np.ones_like(Qpts) + def convert_to_fiat(self, qpts, pts, wts): + pt_dict = {tuple(p) : [(w, tuple())] for (p, w) in zip(pts.T, wts)} + return pt_dict, {} + def __repr__(self): return "H1" @@ -86,6 +90,15 @@ def plot(self, ax, coord, trace_entity, g, **kwargs): vec = np.cross(basis[0], basis[1]) ax.quiver(*coord, *vec, **kwargs) + def convert_to_fiat(self, qpts, pts, wts): + f_at_qpts = pts + shp = tuple(f_at_qpts.shape[:-1]) + weights = np.transpose(np.multiply(f_at_qpts, wts), (-1,) + tuple(range(len(shp)))) + alphas = list(np.ndindex(shp)) + pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)} + return pt_dict, {} + + def tabulate(self, Qpts, trace_entity, g): entityBasis = np.array(trace_entity.basis_vectors()) cellEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) @@ -122,6 +135,14 @@ def tabulate(self, Qpts, trace_entity, g): subEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) result = np.matmul(tangent, subEntityBasis) return result + + def convert_to_fiat(self, qpts, pts, wts): + f_at_qpts = pts + shp = tuple(f_at_qpts.shape[:-1]) + weights = np.transpose(np.multiply(f_at_qpts, wts), (-1,) + tuple(range(len(shp)))) + alphas = list(np.ndindex(shp)) + pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)} + return pt_dict, {} def plot(self, ax, coord, trace_entity, g, **kwargs): permuted = self.domain.permute_entities(g, trace_entity.dimension) diff --git a/fuse/triples.py b/fuse/triples.py index f977b2a..78d5085 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -122,7 +122,19 @@ def to_fiat(self): for i in range(len(dofs)): if entity[1] == dofs[i].trace_entity.id - min_ids[dim]: entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].append(counter) - nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) + if hasattr(dofs[i], "convert_to_fiat_new"): + nodes.append(dofs[i].convert_to_fiat_new(ref_el, degree)) + print("old") + + print(dofs[i].convert_to_fiat(ref_el, degree).pt_dict) + print(dofs[i].convert_to_fiat(ref_el, degree).target_shape) + print("new") + + print(dofs[i].convert_to_fiat_new(ref_el, degree).pt_dict) + print(dofs[i].convert_to_fiat_new(ref_el, degree).target_shape) + else: + raise ValueError("using old") + nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) counter += 1 entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) self.matrices = self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set) @@ -198,7 +210,7 @@ def compute_dense_matrix(self, ref_el, entity_ids, nodes, poly_set): A = dualmat.reshape((shp[0], -1)) B = old_coeffs.reshape((shp[0], -1)) V = np.dot(A, np.transpose(B)) - + print(V) with warnings.catch_warnings(): warnings.filterwarnings("error") try: diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index ee05805..46c929e 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -578,7 +578,7 @@ def test_project_3d(elem_gen, elem_code, deg): assert np.allclose(out.dat.data, f.dat.data, rtol=1e-5) - +@pytest.mark.xfail(reason='Derivative nodes to fiat') def test_create_hermite(): deg = 3 cell = polygon(3) From c6df7bf4146a2d8cbb003652c014cb2e619b6f88 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 28 Jan 2025 17:01:32 +0000 Subject: [PATCH 06/10] convert derivative nodes to fiat def, lint --- fuse/cells.py | 5 +---- fuse/dof.py | 20 ++++++++------------ fuse/traces.py | 29 +++++++++++++++++++++++------ fuse/triples.py | 19 ++++--------------- test/test_convert_to_fiat.py | 23 +++++++++++------------ test/test_dofs.py | 25 ++++++++++++++++++++++++- 6 files changed, 71 insertions(+), 50 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index b487606..9bcfdf1 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -798,8 +798,6 @@ def __init__(self, cell, name=None): name = "FuseCell" self.name = name - - # verts = [cell.get_node(v, return_coords=True) for v in cell.ordered_vertices()] verts = cell.vertices(return_coords=True) topology = cell.get_topology() @@ -942,7 +940,7 @@ def to_fiat(self): return self.cell_complex.to_fiat(name=self.cellname()) def __repr__(self): - return super(CellComplexToUFL, self).__repr__() + return super(CellComplexToUFL, self).__repr__() def reconstruct(self, **kwargs): """Reconstruct this cell, overwriting properties by those in kwargs.""" @@ -980,4 +978,3 @@ def constructCellComplex(name): return TensorProductPoint(*components).to_ufl(name) else: raise TypeError("Cell complex construction undefined for {}".format(str(name))) - diff --git a/fuse/dof.py b/fuse/dof.py index 65e408f..edcd672 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -37,7 +37,7 @@ def __call__(self, kernel, v, cell): def convert_to_fiat(self, ref_el, dof, interpolant_deg): pt = dof.eval(MyTestFunction(lambda *x: x)) return PointEvaluation(ref_el, pt) - + def get_pts(self, ref_el, total_degree): entity = ref_el.construct_subelement(self.entity.dim()) return [(0,) * entity.get_spatial_dimension()], [1], 1 @@ -272,24 +272,20 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non self.sub_id = generator_id def convert_to_fiat(self, ref_el, interpolant_degree): - return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree) - - def convert_to_fiat_new(self, ref_el, interpolant_degree): total_degree = self.kernel.degree() + interpolant_degree pts, wts, jdet = self.pairing.get_pts(ref_el, total_degree) f_pts = self.kernel.tabulate(pts).T / jdet # TODO need to work out how i can discover the shape in a better way if isinstance(self.pairing, DeltaPairing): shp = tuple() - pt_dict = {tuple(p) : [(w, tuple())] for (p, w) in zip(f_pts.T, wts)} + pt_dict = {tuple(p): [(w, tuple())] for (p, w) in zip(f_pts.T, wts)} else: shp = tuple(f_pts.shape[:-1]) weights = np.transpose(np.multiply(f_pts, wts), (-1,) + tuple(range(len(shp)))) alphas = list(np.ndindex(shp)) pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(pts, weights)} - return Functional(ref_el, shp, pt_dict, {}, self.__repr__()) - + return [Functional(ref_el, shp, pt_dict, {}, self.__repr__())] def __repr__(self, fn="v"): return str(self.pairing).format(fn=fn, kernel=self.kernel) @@ -329,16 +325,16 @@ def tabulate(self, Qpts): immersion = self.target_space.tabulate(Qpts, self.trace_entity, self.g) res = self.kernel.tabulate(Qpts) return immersion*res - - def convert_to_fiat_new(self, ref_el, interpolant_degree): + + def convert_to_fiat(self, ref_el, interpolant_degree): total_degree = self.kernel.degree() + interpolant_degree pts, wts, jdet = self.pairing.get_pts(ref_el, total_degree) f_pts = self.kernel.tabulate(pts, self.attachment) attached_pts = [self.attachment(*p) for p in pts] immersion = self.target_space.tabulate(f_pts, self.trace_entity, self.g) - + f_pts = (f_pts * immersion).T / jdet - pt_dict, deriv_dict = self.target_space.convert_to_fiat(attached_pts, f_pts, wts) + dict_list = self.target_space.convert_to_fiat(attached_pts, f_pts, wts) # breakpoint() # TODO need to work out how i can discover the shape in a better way @@ -346,7 +342,7 @@ def convert_to_fiat_new(self, ref_el, interpolant_degree): shp = tuple() else: shp = tuple(f_pts.shape[:-1]) - return Functional(ref_el, shp, pt_dict, deriv_dict, self.__repr__()) + return [Functional(ref_el, shp, pt_dict, deriv_dict, self.__repr__()) for (pt_dict, deriv_dict) in dict_list] def __call__(self, g): permuted = self.cell.permute_entities(g, self.trace_entity.dim()) diff --git a/fuse/traces.py b/fuse/traces.py index 3d9f07b..7aead3b 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -55,8 +55,8 @@ def tabulate(self, Qpts, trace_entity, g): return np.ones_like(Qpts) def convert_to_fiat(self, qpts, pts, wts): - pt_dict = {tuple(p) : [(w, tuple())] for (p, w) in zip(pts.T, wts)} - return pt_dict, {} + pt_dict = {tuple(p): [(w, tuple())] for (p, w) in zip(pts.T, wts)} + return [(pt_dict, {})] def __repr__(self): return "H1" @@ -96,8 +96,7 @@ def convert_to_fiat(self, qpts, pts, wts): weights = np.transpose(np.multiply(f_at_qpts, wts), (-1,) + tuple(range(len(shp)))) alphas = list(np.ndindex(shp)) pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)} - return pt_dict, {} - + return [(pt_dict, {})] def tabulate(self, Qpts, trace_entity, g): entityBasis = np.array(trace_entity.basis_vectors()) @@ -135,14 +134,14 @@ def tabulate(self, Qpts, trace_entity, g): subEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) result = np.matmul(tangent, subEntityBasis) return result - + def convert_to_fiat(self, qpts, pts, wts): f_at_qpts = pts shp = tuple(f_at_qpts.shape[:-1]) weights = np.transpose(np.multiply(f_at_qpts, wts), (-1,) + tuple(range(len(shp)))) alphas = list(np.ndindex(shp)) pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)} - return pt_dict, {} + return [(pt_dict, {})] def plot(self, ax, coord, trace_entity, g, **kwargs): permuted = self.domain.permute_entities(g, trace_entity.dimension) @@ -179,10 +178,28 @@ def apply(*x): return tuple(result) return apply + def convert_to_fiat(self, qpts, pts, wts): + shp = tuple(pts.shape[0:]) + alphas = [] + for i in range(pts.shape[0]): + new = np.zeros_like(shp) + new[i] = 1 + alphas += [tuple(new)] + deriv_dicts = [] + for alpha in alphas: + deriv_dicts += [{tuple(p): [(1.0, tuple(alpha), tuple())] for p in pts.T}] + + # self.alpha = tuple(alpha) + # self.order = sum(self.alpha) + return [({}, d) for d in deriv_dicts] + def plot(self, ax, coord, trace_entity, g, **kwargs): circle1 = plt.Circle(coord, 0.075, fill=False, **kwargs) ax.add_patch(circle1) + def tabulate(self, Qpts, trace_entity, g): + return np.ones_like(Qpts) + def __repr__(self): return "Grad" diff --git a/fuse/triples.py b/fuse/triples.py index 78d5085..5b944c6 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -121,21 +121,10 @@ def to_fiat(self): dim = entity[0] for i in range(len(dofs)): if entity[1] == dofs[i].trace_entity.id - min_ids[dim]: - entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].append(counter) - if hasattr(dofs[i], "convert_to_fiat_new"): - nodes.append(dofs[i].convert_to_fiat_new(ref_el, degree)) - print("old") - - print(dofs[i].convert_to_fiat(ref_el, degree).pt_dict) - print(dofs[i].convert_to_fiat(ref_el, degree).target_shape) - print("new") - - print(dofs[i].convert_to_fiat_new(ref_el, degree).pt_dict) - print(dofs[i].convert_to_fiat_new(ref_el, degree).target_shape) - else: - raise ValueError("using old") - nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) - counter += 1 + fiat_dofs = dofs[i].convert_to_fiat(ref_el, degree) + nodes.extend(fiat_dofs) + entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].extend([counter + i for i in range(len(fiat_dofs))]) + counter += len(fiat_dofs) entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) self.matrices = self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set) form_degree = 1 if self.spaces[0].set_shape else 0 diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 46c929e..b89317f 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -154,6 +154,7 @@ def create_cg2_tri(cell): DOFGenerator(edge_xs, C3, S1)]) return cg + def create_hermite(tri): vert = tri.vertices()[0] @@ -166,14 +167,11 @@ def create_hermite(tri): v_derv_xs = [immerse(tri, dg0, TrGrad)] v_derv_dofs = DOFGenerator(v_derv_xs, S3/S2, S1) - v_derv2_xs = [immerse(tri, dg0, TrHess)] - v_derv2_dofs = DOFGenerator(v_derv2_xs, S3/S2, S1) - i_xs = [DOF(DeltaPairing(), PointKernel((0, 0)))] i_dofs = DOFGenerator(i_xs, S1, S1) her = ElementTriple(tri, (P3, CellH2, C0), - [v_dofs, v_derv_dofs, v_derv2_dofs, i_dofs]) + [v_dofs, v_derv_dofs, i_dofs]) return her @@ -578,7 +576,8 @@ def test_project_3d(elem_gen, elem_code, deg): assert np.allclose(out.dat.data, f.dat.data, rtol=1e-5) -@pytest.mark.xfail(reason='Derivative nodes to fiat') + +@pytest.mark.xfail(reason='Handling generation of multiple fiat nodes from one in permutations') def test_create_hermite(): deg = 3 cell = polygon(3) @@ -595,16 +594,16 @@ def test_create_hermite(): Qpts, _ = Q.get_points(), Q.get_weights() fiat_vals = fiat_elem.tabulate(0, Qpts) - # my_vals = my_elem.tabulate(0, Qpts) + my_vals = my_elem.tabulate(0, Qpts) fiat_vals = flatten(fiat_vals[(0,) * sd]) - # my_vals = flatten(my_vals[(0,) * sd]) + my_vals = flatten(my_vals[(0,) * sd]) - # (x, res, _, _) = np.linalg.lstsq(fiat_vals.T, my_vals.T) - # x1 = np.linalg.inv(x) - # assert np.allclose(np.linalg.norm(my_vals.T - fiat_vals.T @ x), 0) - # assert np.allclose(np.linalg.norm(fiat_vals.T - my_vals.T @ x1), 0) - # assert np.allclose(res, 0) + (x, res, _, _) = np.linalg.lstsq(fiat_vals.T, my_vals.T) + x1 = np.linalg.inv(x) + assert np.allclose(np.linalg.norm(my_vals.T - fiat_vals.T @ x), 0) + assert np.allclose(np.linalg.norm(fiat_vals.T - my_vals.T @ x1), 0) + assert np.allclose(res, 0) def test_investigate_dpc(): diff --git a/test/test_dofs.py b/test/test_dofs.py index eb1c046..687c82c 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -1,5 +1,5 @@ from fuse import * -from test_convert_to_fiat import create_cg1, create_dg1, construct_cg3, construct_rt, construct_nd +from test_convert_to_fiat import create_cg1, create_dg1, construct_cg3, construct_rt, construct_nd, create_hermite import sympy as sp import numpy as np @@ -106,3 +106,26 @@ def test_permute_nd(): # print(g) # print(nd.cell.permute_entities(g, 0)) # print(nd.cell.permute_entities(g, 1)) + + +def test_convert_dofs(): + cell = polygon(3) + + cg3 = create_hermite(cell) + + for dof in cg3.generate(): + print(dof) + # print("old") + # # old = dof.convert_to_fiat(cell.to_fiat(), 5).pt_dict + # print(old) + # print("new") + new = dof.convert_to_fiat_new(cell.to_fiat(), 5)[0].pt_dict + print(new) + new = dof.convert_to_fiat_new(cell.to_fiat(), 5)[0].deriv_dict + print(new) + from FIAT.hermite import CubicHermite + fiat_elem = CubicHermite(cell.to_fiat(), 3) + + print(len([n.pt_dict for n in fiat_elem.dual.nodes])) + print([n.pt_dict for n in fiat_elem.dual.nodes]) + print([n.deriv_dict for n in fiat_elem.dual.nodes]) From a508380969b4fe5ad923ac1c5041cc23e3a3b350 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 18 Mar 2025 09:52:36 +0000 Subject: [PATCH 07/10] fixing post rebase, draft of bfs element --- fuse/cells.py | 4 ---- fuse/triples.py | 9 ++++++--- test/bfs.py | 23 +++++++++++++++++++++++ 3 files changed, 29 insertions(+), 7 deletions(-) create mode 100644 test/bfs.py diff --git a/fuse/cells.py b/fuse/cells.py index 9bcfdf1..c188de0 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -383,17 +383,13 @@ def get_topology(self): self.topology = {} self.topology_unrelabelled = {} - self.topology_unrelabelled = {} for i in range(len(structure)): dimension = structure[i] self.topology[i] = {} self.topology_unrelabelled[i] = {} - self.topology_unrelabelled[i] = {} for node in dimension: self.topology[i][node - min_ids[i]] = tuple([relabelled_verts[vert] for vert in self.get_node(node).ordered_vertices()]) self.topology_unrelabelled[i][node - min_ids[i]] = tuple([vert - min_ids[0] for vert in self.get_node(node).ordered_vertices()]) - return self.topology_unrelabelled - self.topology_unrelabelled[i][node - min_ids[i]] = tuple([vert - min_ids[0] for vert in self.get_node(node).ordered_vertices()]) return self.topology_unrelabelled def get_starter_ids(self): diff --git a/fuse/triples.py b/fuse/triples.py index 5b944c6..fee2079 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -125,18 +125,21 @@ def to_fiat(self): nodes.extend(fiat_dofs) entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].extend([counter + i for i in range(len(fiat_dofs))]) counter += len(fiat_dofs) + print("nodes", nodes) entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) - self.matrices = self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set) + # self.matrices = self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set) form_degree = 1 if self.spaces[0].set_shape else 0 print("my", [n.pt_dict for n in nodes]) print(entity_perms) print(entity_ids) - print(ref_el.vertices) + print(nodes) print() # TODO: Change this when Dense case in Firedrake if pure_perm: + self.matrices = self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set) dual = DualSet(nodes, ref_el, entity_ids, entity_perms) else: + self.matrices = None dual = DualSet(nodes, ref_el, entity_ids) return CiarletElement(poly_set, dual, degree, form_degree) @@ -221,7 +224,7 @@ def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): if g.perm.is_Identity: res_dict[dim][e_id][val] = np.eye(len(nodes)) else: - new_nodes = [d(g).convert_to_fiat(ref_el, degree) for d in self.generate()] + new_nodes = sum([d(g).convert_to_fiat(ref_el, degree) for d in self.generate()],[]) transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_ids, new_nodes, poly_set) res_dict[dim][e_id][val] = np.matmul(transformed_basis, original_V.T) return res_dict diff --git a/test/bfs.py b/test/bfs.py new file mode 100644 index 0000000..f09344b --- /dev/null +++ b/test/bfs.py @@ -0,0 +1,23 @@ +from firedrake import * +from fuse import * +import numpy as np +from test_convert_to_fiat import create_dg1 +from test_tensor_prod import mass_solve + +def her_int(): + deg = 3 + cell = Point(1, [Point(0), Point(0)], vertex_num=2) + vert_dg = create_dg1(cell.vertices()[0]) + xs = [immerse(cell, vert_dg, TrH1), immerse(cell, vert_dg, TrGrad)] + + Pk = PolynomialSpace(deg) + her = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, S2, S1)) + return her + +r = 3 +her1 = her_int() +her2 = her_int() +bfs = tensor_product(her1, her2).flatten() +mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) +U = FunctionSpace(mesh, bfs.to_ufl()) +mass_solve(U) \ No newline at end of file From 7da2601e7dafaf4afc9230637c6f19367d789694 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Sat, 22 Mar 2025 14:43:15 +0000 Subject: [PATCH 08/10] lint --- fuse/triples.py | 2 +- test/bfs.py | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/fuse/triples.py b/fuse/triples.py index fee2079..c778737 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -224,7 +224,7 @@ def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): if g.perm.is_Identity: res_dict[dim][e_id][val] = np.eye(len(nodes)) else: - new_nodes = sum([d(g).convert_to_fiat(ref_el, degree) for d in self.generate()],[]) + new_nodes = sum([d(g).convert_to_fiat(ref_el, degree) for d in self.generate()], []) transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_ids, new_nodes, poly_set) res_dict[dim][e_id][val] = np.matmul(transformed_basis, original_V.T) return res_dict diff --git a/test/bfs.py b/test/bfs.py index f09344b..e3e2e24 100644 --- a/test/bfs.py +++ b/test/bfs.py @@ -1,9 +1,9 @@ from firedrake import * from fuse import * -import numpy as np from test_convert_to_fiat import create_dg1 from test_tensor_prod import mass_solve + def her_int(): deg = 3 cell = Point(1, [Point(0), Point(0)], vertex_num=2) @@ -14,10 +14,11 @@ def her_int(): her = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, S2, S1)) return her + r = 3 her1 = her_int() her2 = her_int() bfs = tensor_product(her1, her2).flatten() mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) U = FunctionSpace(mesh, bfs.to_ufl()) -mass_solve(U) \ No newline at end of file +mass_solve(U) From 5dca4dfabc31751dec7779abea86f45d058c66fc Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 24 Apr 2025 11:47:05 +0100 Subject: [PATCH 09/10] Derivative type nodes seem functional --- fuse/traces.py | 4 ++-- fuse/triples.py | 8 ++++++-- test/bfs.py | 3 ++- test/test_convert_to_fiat.py | 6 ++++-- test/test_dofs.py | 31 +++++++++++++++++++------------ 5 files changed, 33 insertions(+), 19 deletions(-) diff --git a/fuse/traces.py b/fuse/traces.py index 6da204a..81fa713 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -180,10 +180,10 @@ def apply(*x): return apply def convert_to_fiat(self, qpts, pts, wts): - shp = tuple(pts.shape[0:]) + shp = (self.domain.get_spatial_dimension(),) alphas = [] for i in range(pts.shape[0]): - new = np.zeros_like(shp) + new = np.zeros(shp, dtype=int) new[i] = 1 alphas += [tuple(new)] deriv_dicts = [] diff --git a/fuse/triples.py b/fuse/triples.py index 1509905..498159d 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -128,10 +128,12 @@ def to_fiat(self): entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) # self.matrices = self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set) form_degree = 1 if self.spaces[0].set_shape else 0 - print("my", [n.pt_dict for n in nodes]) + print("pts", [n.pt_dict for n in nodes]) + print("derivs", [n.deriv_dict for n in nodes]) print(entity_perms) print(entity_ids) print(nodes) + print(len(nodes)) print() # TODO: Change this when Dense case in Firedrake if pure_perm: @@ -140,6 +142,7 @@ def to_fiat(self): else: self.matrices = None dual = DualSet(nodes, ref_el, entity_ids) + # breakpoint() return CiarletElement(poly_set, dual, degree, form_degree) def to_tikz(self, show=True, scale=3): @@ -290,7 +293,8 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): if pure_perm is False: # TODO think about where this call goes - return self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set), pure_perm + return None, False + # return self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set), pure_perm dof_id_mat = np.eye(len(dofs)) oriented_mats_by_entity = {} diff --git a/test/bfs.py b/test/bfs.py index e3e2e24..a03c3b2 100644 --- a/test/bfs.py +++ b/test/bfs.py @@ -11,7 +11,8 @@ def her_int(): xs = [immerse(cell, vert_dg, TrH1), immerse(cell, vert_dg, TrGrad)] Pk = PolynomialSpace(deg) - her = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, S2, S1)) + # TODO fix faking out permutations + her = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, S2, S2)) return her diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index b89317f..2cd112d 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -165,7 +165,8 @@ def create_hermite(tri): v_dofs = DOFGenerator(v_xs, S3/S2, S1) v_derv_xs = [immerse(tri, dg0, TrGrad)] - v_derv_dofs = DOFGenerator(v_derv_xs, S3/S2, S1) + # TODO fix perms so this can be S1 + v_derv_dofs = DOFGenerator(v_derv_xs, S3/S2, S3) i_xs = [DOF(DeltaPairing(), PointKernel((0, 0)))] i_dofs = DOFGenerator(i_xs, S1, S1) @@ -505,6 +506,7 @@ def test_non_tensor_quad(): (lambda cell: CR_n(cell, 3), "CR", 1), (create_cf, "CR", 1), # Don't think Crouzeix Falk in in Firedrake (construct_cg3, "CG", 3), + (create_hermite, "HER", 3), pytest.param(construct_nd, "N1curl", 1, marks=pytest.mark.xfail(reason='Dense Matrices needed')), pytest.param(construct_rt, "RT", 1, marks=pytest.mark.xfail(reason='Dense Matrices needed'))]) def test_project(elem_gen, elem_code, deg): @@ -577,7 +579,7 @@ def test_project_3d(elem_gen, elem_code, deg): assert np.allclose(out.dat.data, f.dat.data, rtol=1e-5) -@pytest.mark.xfail(reason='Handling generation of multiple fiat nodes from one in permutations') +# @pytest.mark.xfail(reason='Handling generation of multiple fiat nodes from one in permutations') def test_create_hermite(): deg = 3 cell = polygon(3) diff --git a/test/test_dofs.py b/test/test_dofs.py index 96bbb57..4028d79 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -113,19 +113,26 @@ def test_convert_dofs(): cg3 = create_hermite(cell) - for dof in cg3.generate(): - print(dof) - # print("old") - # # old = dof.convert_to_fiat(cell.to_fiat(), 5).pt_dict - # print(old) - # print("new") - new = dof.convert_to_fiat_new(cell.to_fiat(), 5)[0].pt_dict - print(new) - new = dof.convert_to_fiat_new(cell.to_fiat(), 5)[0].deriv_dict - print(new) + # for dof in cg3.generate(): + # print(dof) + # # print("old") + # # # old = dof.convert_to_fiat(cell.to_fiat(), 5).pt_dict + # # print(old) + # # print("new") + # new = dof.convert_to_fiat(cell.to_fiat(), 5)[0].pt_dict + # print(new) + # new = dof.convert_to_fiat(cell.to_fiat(), 5)[0].deriv_dict + # print(new) + # print(len(cg3.generate())) + fuse_elem = cg3.to_fiat() + for n in fuse_elem.dual.nodes: + print(n.pt_dict) + print(n.deriv_dict) + from FIAT.hermite import CubicHermite fiat_elem = CubicHermite(cell.to_fiat(), 3) print(len([n.pt_dict for n in fiat_elem.dual.nodes])) - print([n.pt_dict for n in fiat_elem.dual.nodes]) - print([n.deriv_dict for n in fiat_elem.dual.nodes]) + for n in fiat_elem.dual.nodes: + print(n.pt_dict) + print(n.deriv_dict) From 78ac5be8d5dc1b01af63b566c9bd9edf6a04c261 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 1 May 2025 13:34:11 +0100 Subject: [PATCH 10/10] mass solve bfs in, not working --- fuse/triples.py | 1 - test/bfs.py | 2 +- test/test_dofs.py | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/fuse/triples.py b/fuse/triples.py index 498159d..2de1bcb 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -142,7 +142,6 @@ def to_fiat(self): else: self.matrices = None dual = DualSet(nodes, ref_el, entity_ids) - # breakpoint() return CiarletElement(poly_set, dual, degree, form_degree) def to_tikz(self, show=True, scale=3): diff --git a/test/bfs.py b/test/bfs.py index a03c3b2..1152ec5 100644 --- a/test/bfs.py +++ b/test/bfs.py @@ -22,4 +22,4 @@ def her_int(): bfs = tensor_product(her1, her2).flatten() mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) U = FunctionSpace(mesh, bfs.to_ufl()) -mass_solve(U) +mass_solve(U) \ No newline at end of file diff --git a/test/test_dofs.py b/test/test_dofs.py index 4028d79..3b964a7 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -135,4 +135,4 @@ def test_convert_dofs(): print(len([n.pt_dict for n in fiat_elem.dual.nodes])) for n in fiat_elem.dual.nodes: print(n.pt_dict) - print(n.deriv_dict) + print(n.deriv_dict) \ No newline at end of file