From a1c7558cc771442da7affe56c0725579bf105338 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 17 Oct 2025 17:39:58 +0100 Subject: [PATCH 01/12] issues with poly degree and adding parameterisation kernels --- fuse/cells.py | 27 +++++ fuse/dof.py | 69 ++++++++++-- fuse/groups.py | 6 +- fuse/spaces/polynomial_spaces.py | 6 +- fuse/triples.py | 9 +- test/test_dofs.py | 35 ++++++ test/test_orientations.py | 178 ++++++++++++++++++++----------- 7 files changed, 254 insertions(+), 76 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index fb9694a..b9b03db 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -656,6 +656,29 @@ def to_tikz(self, show=True, scale=3): return "\n".join(tikz_commands) return tikz_commands + + def generate_facet_parameterisation(self, facet_num): + #facet = self.d_entities(self.dimension - 1)[facet_num] + facet = self.get_node(facet_num) + facet_dim = facet.dimension + if facet_dim != self.dimension - 1: + raise ValueError(f"Supplied node {facet_num} is not a facet") + if facet_dim > 1: + raise NotImplementedError("Facet parameterisation is not implemented for dimensions greater than 1") + verts = facet.vertices() + v_coords = np.array([self.get_node(v.id, return_coords=True) for v in verts]) + midpoint = np.average(v_coords, axis=0) + stacked = np.c_[np.ones((self.dimension,)), v_coords[:, 0].reshape(self.dimension,1)] + b = np.array([0, 1]) + coeffs = np.linalg.solve(stacked, b) + + symbol_names = ["x", "y", "z"] + symbols = [1] + [sp.Symbol(symbol_names[d]) for d in range(facet_dim)] + res = 0 + for d in range(facet_dim + 1): + res += coeffs[d] * symbols[d] + return res + def plot(self, show=True, plain=False, ax=None, filename=None): """ for now into 2 dimensional space """ if self.dimension == 3: @@ -766,6 +789,10 @@ def attachment(self, source, dst): assert all(np.isclose(val, vals[0]).all() for val in vals) return lambda *x: fold_reduce(attachments[0], *x) + + + + def cell_attachment(self, dst): if not isinstance(dst, int): diff --git a/fuse/dof.py b/fuse/dof.py index 52afcb5..e51ccdf 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -152,8 +152,13 @@ def _from_dict(obj_dict): class BaseKernel(): def __init__(self): + self.cell = None + self.entity = None self.attachment = False + def add_context(self, cell, entity): + return self + def permute(self, g): raise NotImplementedError("This method should be implemented by the subclass") @@ -219,9 +224,9 @@ def degree(self, interpolant_degree): return self.fn.as_poly().total_degree() * interpolant_degree def permute(self, g): - return self - # new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) - # return PolynomialKernel(new_fn, symbols=self.syms) + return self + #new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) + #return PolynomialKernel(new_fn, symbols=self.syms) def __call__(self, *args): res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) @@ -230,9 +235,6 @@ def __call__(self, *args): return res def tabulate(self, Qpts, attachment=None): - # TODO do we need to attach qpts - # if attachment: - # return np.array([self(*attachment(*pt)) for pt in Qpts]).astype(np.float64) return np.array([self(*pt) for pt in Qpts]).astype(np.float64) def _to_dict(self): @@ -277,7 +279,58 @@ def dict_id(self): def _from_dict(obj_dict): return ComponentKernel(obj_dict["comp"]) +class ParameterisationKernel(BaseKernel): + + def __init__(self, g=None): + self.g = g + self.fn = None + super(ParameterisationKernel, self).__init__() + + def add_context(self, cell, entity, fn=None): + new_cls = ParameterisationKernel() + new_cls.cell = cell + new_cls.entity = entity + new_cls.fn = cell.generate_facet_parameterisation(entity.id) + return new_cls + + def __repr__(self): + return f"[{self.comp}]" + + def degree(self, interpolant_degree): + return interpolant_degree + + def permute(self, g): + new_cls = ParameterisationKernel(g) + if self.cell is not None: + new_cls.add_context(self.cell, self.entity) + return new_cls + + def __call__(self, *args): + if self.fn is None: + raise ArgumentError("Function not defined") + + if self.g is not None: + self.fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) + res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) + if not hasattr(res, '__iter__'): + return [res] + return res + + def tabulate(self, Qpts, attachment=None): + + return np.array([1 for _ in Qpts]).astype(np.float64) + return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + + def _to_dict(self): + o_dict = {"comp": self.comp} + return o_dict + + def dict_id(self): + return "ParameterisationKernel" + + def _from_dict(obj_dict): + return ParameterisationKernel(obj_dict["comp"]) class DOF(): @@ -299,6 +352,8 @@ def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=N self.generation = generation if entity is not None: self.pairing = self.pairing.add_entity(entity) + if self.trace_entity is not None and self.cell is not None: + self.kernel = kernel.add_context(cell, entity) def __call__(self, g): new_generation = self.generation.copy() @@ -430,3 +485,5 @@ def _to_dict(self): def dict_id(self): return "Function" + + diff --git a/fuse/groups.py b/fuse/groups.py index a08d96e..14a5a41 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -9,12 +9,12 @@ def perm_matrix_to_perm_array(p_mat): summed = np.sum(p_mat, axis=0) - if np.all(summed == np.zeros_like(summed)): + if np.allclose(summed, np.zeros_like(summed)): return list(np.zeros_like(summed)) - assert np.all(summed == np.ones_like(summed)) + assert np.allclose(summed, np.ones_like(summed)) res = [] for row in p_mat: - indices = list(row).index(1) + indices = list(np.isclose(row, 1)).index(True) res += [indices] return res diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index c05844a..56bbce1 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -147,9 +147,11 @@ 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] - maxdegree = max([space.maxdegree for space in spaces]) - mindegree = min([space.mindegree for space in spaces]) + 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)]) vec = any([s.set_shape for s in spaces]) super(ConstructedPolynomialSpace, self).__init__(maxdegree, -1, mindegree, set_shape=vec) diff --git a/fuse/triples.py b/fuse/triples.py index 3d4e3f7..63206bc 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -353,9 +353,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): dofs = self.generate() min_ids = self.cell.get_starter_ids() entity_associations, pure_perm, sub_pure_perm = self._entity_associations(dofs) - if pure_perm is False: + #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), None, pure_perm + #return self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set), None, pure_perm + # return self.matrices_by_entity, None, pure_perm oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(dofs) # for each entity, look up generation on that entity and permute the @@ -436,8 +437,8 @@ def orient_mat_perms(self): ents = self.cell.d_entities(dim) for e in ents: e_id = e.id - min_ids[dim] - members = e.group.members() - if entity_orientations[num_ents + e_id] != 0 and dim < self.cell.dim(): + members = e.group.members()#and dim < self.cell.dim() + if entity_orientations[num_ents + e_id] != 0 : modifier = self.matrices[dim][e_id][entity_orientations[num_ents+e_id]] reverse_modifier_val = (~e.group.get_member_by_val(entity_orientations[num_ents+e_id])).numeric_rep() reverse_modifier = self.matrices[dim][e_id][reverse_modifier_val] diff --git a/test/test_dofs.py b/test/test_dofs.py index 5b56f43..ed81804 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -138,3 +138,38 @@ def test_permute_nodes(): print([n.pt_dict for n in nodes]) for g in cg1.cell.group.members(): print(g, [d(g).convert_to_fiat(cell.to_fiat(), degree).pt_dict for d in cg1.generate()]) + +def test_edge_parametrisation(): + tri = polygon(3) + for i in tri.d_entities(1): + print(tri.generate_facet_parameterisation(i.id)) + from fuse.dof import ParameterisationKernel + + + deg = 2 + edge = tri.edges()[0] + x = sp.Symbol("x") + y = sp.Symbol("y") + + xs = [DOF(L2Pairing(), ParameterisationKernel())] + + + dofs = DOFGenerator(xs, S2, S2) + int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) + + + xs = [DOF(L2Pairing(), ComponentKernel((0,))), + DOF(L2Pairing(), ComponentKernel((1,)))] + center_dofs = DOFGenerator(xs, S1, S3) + xs = [immerse(tri, int_ned1, TrHCurl)] + tri_dofs = DOFGenerator(xs, C3, S1) + + vec_Pk = PolynomialSpace(deg - 1, set_shape=True) + Pk = PolynomialSpace(deg - 1) + M = sp.Matrix([[y, -x]]) + nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M + + + ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) + for n in ned.generate(): + print(n) diff --git a/test/test_orientations.py b/test/test_orientations.py index de81546..026c7ec 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -1,10 +1,12 @@ import unittest.mock as mock +import pytest from firedrake import * from fuse import * import sympy as sp -from test_convert_to_fiat import create_cg1, helmholtz_solve, construct_nd, construct_rt +from test_convert_to_fiat import create_cg1, helmholtz_solve, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg1 import os +#with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): def dummy_dof_perms(cls, *args, **kwargs): # return -1s of right shape here oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) @@ -19,18 +21,96 @@ def dummy_dof_perms(cls, *args, **kwargs): oriented_mats_by_entity[key1][key2][key3] = 100 * np.identity(val3.shape[0]) return oriented_mats_by_entity, False, None +def construct_nd2(tri=None): + if tri is None: + tri = polygon(3) + deg = 2 + edge = tri.edges()[0] + x = sp.Symbol("x") + y = sp.Symbol("y") + + #xs = [DOF(L2Pairing(), PolynomialKernel((x + 1)/2, symbols=[x])), + # DOF(L2Pairing(), PolynomialKernel(1))] + xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), + DOF(L2Pairing(), PolynomialKernel((-np.sqrt(2)/np.sqrt(3))*(6*x + 3), symbols=[x]))] + + + dofs = DOFGenerator(xs, S1, S2) + int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) + + + xs = [DOF(L2Pairing(), ComponentKernel((0,))), + DOF(L2Pairing(), ComponentKernel((1,)))] + center_dofs = DOFGenerator(xs, S1, S3) + xs = [immerse(tri, int_ned1, TrHCurl)] + tri_dofs = DOFGenerator(xs, C3, S1) + + vec_Pk = PolynomialSpace(deg - 1, set_shape=True) + Pk = PolynomialSpace(deg - 1) + M = sp.Matrix([[y, -x]]) + nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M + print(nd.degree()) + + + ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) + return ned + +def construct_nd2_for_fiat(tri=None): + if tri is None: + tri = polygon(3) + deg = 2 + edge = tri.edges()[0] + x = sp.Symbol("x") + y = sp.Symbol("y") + + xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), + DOF(L2Pairing(), PolynomialKernel((3*x*np.sqrt(2)/np.sqrt(3)), symbols=[x]))] + + + dofs = DOFGenerator(xs, S1, S2) + int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) + + xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), + DOF(L2Pairing(), PolynomialKernel((-np.sqrt(2)/np.sqrt(3))*(6*x + 3), symbols=[x]))] + dofs = DOFGenerator(xs, S1, S2) + int_ned2 = ElementTriple(tri.edges()[0], (P1, CellHCurl, C0), dofs) + + xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), + DOF(L2Pairing(), PolynomialKernel((-np.sqrt(2)/np.sqrt(3))*(6*x - 3), symbols=[x]))] + dofs = DOFGenerator(xs, S1, S2) + int_ned3 = ElementTriple(tri.edges()[0], (P1, CellHCurl, C0), dofs) + + xs = [DOF(L2Pairing(), ComponentKernel((0,))), + DOF(L2Pairing(), ComponentKernel((1,)))] + center_dofs = DOFGenerator(xs, S1, S3) + xs = [immerse(tri, int_ned2, TrHCurl, node=0)] + xs += [immerse(tri, int_ned1, TrHCurl, node=1)] + xs += [immerse(tri, int_ned3, TrHCurl, node=2)] + tri_dofs = DOFGenerator(xs, S1, S1) + + vec_Pk = PolynomialSpace(deg - 1, set_shape=True) + Pk = PolynomialSpace(deg - 1) + M = sp.Matrix([[y, -x]]) + nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M + -def test_interpolation_values(): + ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) + return ned + +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), + (construct_nd2_for_fiat, "N1curl", 2), + (construct_nd2, "N1curl", 2)]) +def test_interpolation_values(elem_gen, elem_code,deg): cell = polygon(3) - elem = construct_nd(cell) + elem = elem_gen(cell) print() - #with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): - mesh = UnitSquareMesh(3, 3) + #mesh = UnitSquareMesh(1, 1) + mesh = UnitTriangleMesh(0) ones = as_vector((0,1)) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: - V = FunctionSpace(mesh, "N1curl", 1) + V = FunctionSpace(mesh, elem_code, deg) print(V.cell_node_list) u = TestFunction(V) @@ -44,7 +124,7 @@ def test_surface_const_nd(): elem = construct_nd(cell) ones = as_vector((0,1)) - for n in range(1, 6): + for n in range(3, 4): mesh = UnitSquareMesh(n, n) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): @@ -144,54 +224,39 @@ def test_interpolate_vs_project(V): print(expect.dat.data) assert np.allclose(f.dat.data, expect.dat.data, atol=1e-06) -def construct_nd2(tri=None): - if tri is None: - tri = polygon(3) - deg = 2 - edge = tri.edges()[0] - x = sp.Symbol("x") - y = sp.Symbol("y") - - #xs = [DOF(L2Pairing(), PointKernel((-np.sqrt(2),)))] - xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x+1), symbols=[x])), - DOF(L2Pairing(), PolynomialKernel((1/2)*(1-x), symbols=[x]))] - - - dofs = DOFGenerator(xs, S1, S2) - int_ned = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - - xs = [DOF(L2Pairing(), ComponentKernel((0,))), - DOF(L2Pairing(), ComponentKernel((1,)))] - center_dofs = DOFGenerator(xs, S1, S3) - xs = [immerse(tri, int_ned, TrHCurl)] - tri_dofs = DOFGenerator(xs, C3, S1) - - vec_Pk = PolynomialSpace(deg - 1, set_shape=True) - Pk = PolynomialSpace(deg - 1) - M = sp.Matrix([[y, -x]]) - nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M - - - ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) - return ned - def test_degree2_interpolation(): cell = polygon(3) - #elem = construct_nd2(cell) + elem = construct_nd2_for_fiat(cell) mesh = UnitSquareMesh(1,1) - #print("fuse") - #V = FunctionSpace(mesh, elem.to_ufl()) - #test_interpolate_vs_project(V) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, "N1curl", 2) + test_interpolate_vs_project(V) + +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), + (create_cg1, "CG", 1), + (create_dg1, "DG", 1), + (construct_cg3, "CG", 3), + ]) +def test_interpolation(elem_gen, elem_code, deg): + cell = polygon(3) + elem = elem_gen(cell) + mesh = UnitSquareMesh(1,1) - print("firedrake") - V = FunctionSpace(mesh, "N1curl", 2) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) test_interpolate_vs_project(V) def test_create_fiat_nd(): cell = polygon(3) nd = construct_nd2(cell) + breakpoint() + nd_fv = construct_nd2_for_fiat(cell) ref_el = cell.to_fiat() sd = ref_el.get_spatial_dimension() deg = 2 @@ -201,27 +266,18 @@ def test_create_fiat_nd(): print("fiat") for f in fiat_elem.dual_basis(): print(f.pt_dict) - + + print("fiat: fuse's version") + for d in nd_fv.generate(): + print(d.convert_to_fiat(ref_el, deg).pt_dict) + print("fuse") for d in nd.generate(): print(d.convert_to_fiat(ref_el, deg).pt_dict) - print(d) - breakpoint() + #print(d) my_elem = nd.to_fiat() + my_elem = nd_fv.to_fiat() + breakpoint() - 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) From 125b0a695524ac2140ce18c4c1322e43625afbe4 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 24 Oct 2025 10:55:32 +0100 Subject: [PATCH 02/12] debugging --- fuse/triples.py | 4 +--- test/test_orientations.py | 35 ++++++++++++++++++++--------------- 2 files changed, 21 insertions(+), 18 deletions(-) diff --git a/fuse/triples.py b/fuse/triples.py index 63206bc..309c399 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -108,8 +108,6 @@ def to_fiat(self): top = ref_el.get_topology() min_ids = self.cell.get_starter_ids() poly_set = self.spaces[0].to_ON_polynomial_set(ref_el) - #from FIAT.nedelec import Nedelec - #poly_set = Nedelec(ref_el, 2).poly_set for dim in sorted(top): entity_ids[dim] = {i: [] for i in top[dim]} @@ -129,7 +127,7 @@ def to_fiat(self): if not pure_perm: self.matrices = mat_perms self.reverse_dof_perms() - self.orient_mat_perms() + #self.orient_mat_perms() form_degree = 1 if self.spaces[0].set_shape else 0 # TODO: Change this when Dense case in Firedrake diff --git a/test/test_orientations.py b/test/test_orientations.py index 026c7ec..5997047 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -3,7 +3,7 @@ from firedrake import * from fuse import * import sympy as sp -from test_convert_to_fiat import create_cg1, helmholtz_solve, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg1 +from test_convert_to_fiat import create_cg1, helmholtz_solve, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg0, create_dg1 import os #with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): @@ -104,14 +104,15 @@ def test_interpolation_values(elem_gen, elem_code,deg): cell = polygon(3) elem = elem_gen(cell) print() - #mesh = UnitSquareMesh(1, 1) - mesh = UnitTriangleMesh(0) + mesh = UnitSquareMesh(1, 1) + #mesh = UnitTriangleMesh(0) ones = as_vector((0,1)) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, elem_code, deg) + print(mesh.entity_orientations) print(V.cell_node_list) u = TestFunction(V) res1= assemble(interpolate(ones, V)) @@ -124,7 +125,7 @@ def test_surface_const_nd(): elem = construct_nd(cell) ones = as_vector((0,1)) - for n in range(3, 4): + for n in range(2, 3): mesh = UnitSquareMesh(n, n) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): @@ -164,33 +165,37 @@ def test_surface_vec(): rt_elem = construct_rt(cell) nd_elem = construct_nd(cell) - for n in range(1, 6): + for n in range(2, 3): mesh = UnitSquareMesh(n, n) x, y = SpatialCoordinate(mesh) normal = FacetNormal(mesh) test_vec = as_vector((-y,x)) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V = FunctionSpace(mesh, rt_elem.to_ufl()) - vec1 = interpolate(test_vec, V) - res1 = assemble(dot(vec1, normal) * ds) - else: - V = FunctionSpace(mesh, "RT", 1) - vec2 = interpolate(test_vec, V) - res1 = assemble(dot(vec2, normal) * ds) - print(f"div {n}: {res1}") + #if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + # V = FunctionSpace(mesh, rt_elem.to_ufl()) + # vec1 = interpolate(test_vec, V) + # res1 = assemble(dot(vec1, normal) * ds) + #else: + # V = FunctionSpace(mesh, "RT", 1) + # vec2 = interpolate(test_vec, V) + # res1 = assemble(dot(vec2, normal) * ds) + #print(f"div {n}: {res1}") if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, nd_elem.to_ufl()) + print(V.cell_node_list) vec1 = interpolate(test_vec, V) + #V2 = FunctionSpace(mesh, create_dg0(cell).to_ufl()) + #u = TestFunction(V2) res2 = assemble(dot(vec1, normal) * ds) else: V = FunctionSpace(mesh, "N1curl", 1) + print(V.cell_node_list) vec1 = interpolate(test_vec, V) res2 = assemble(dot(vec1, normal) * ds) print(f"curl {n}: {res2}") - assert np.allclose(0, res1) + #assert np.allclose(0, res1) assert np.allclose(0, res2) From 3bc51cc0bc9ab9cff8362560d56f632dbb2c2477 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 27 Oct 2025 15:44:04 +0000 Subject: [PATCH 03/12] light refactor --- test/test_orientations.py | 78 +++++++++++++++++++++++---------------- 1 file changed, 46 insertions(+), 32 deletions(-) diff --git a/test/test_orientations.py b/test/test_orientations.py index 5997047..177aff4 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -97,41 +97,55 @@ def construct_nd2_for_fiat(tri=None): ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) return ned -@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), - (construct_nd2_for_fiat, "N1curl", 2), - (construct_nd2, "N1curl", 2)]) +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1)]) +# (construct_nd2_for_fiat, "N1curl", 2), +# (construct_nd2, "N1curl", 2)]) def test_interpolation_values(elem_gen, elem_code,deg): + # this test may not actually make sense proceed with caution + #os.environ["FIREDRAKE_USE_FUSE"] = "1" cell = polygon(3) elem = elem_gen(cell) print() - mesh = UnitSquareMesh(1, 1) - #mesh = UnitTriangleMesh(0) - ones = as_vector((0,1)) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V = FunctionSpace(mesh, elem.to_ufl()) - else: - V = FunctionSpace(mesh, elem_code, deg) - - print(mesh.entity_orientations) - print(V.cell_node_list) - u = TestFunction(V) - res1= assemble(interpolate(ones, V)) - for i in range(len(res1.dat.data)): - print(f"{i}: {res1.dat.data[i]}") - - -def test_surface_const_nd(): + meshes = [UnitTriangleMesh(0), UnitSquareMesh(1,1), UnitSquareMesh(2,2), UnitSquareMesh(3,3)] + vec1 = as_vector((0, 1)) + vec1_expected = [[0, 1, 1], [0, -1, -1, -1, 0], [1, -1, 0,0,1,0,1,1,0,-1,-1,1,0,-1,0,1], [1,-1,0,0,1,0,1,1,0,1,1,-1,0,0,-1,-1,1,1,-1,-1,0,0,-1,-1,0,1,-1,0,1,0,-1,0,1]] + vec2 = as_vector((1, 0)) + vec2_expected = [[1, -1, 0], [1, 0, 1, 0, -1], [0,1,1,1,0,1,0,-1,1,0,1,0,1,1,1,0], [0,1,1,1,0,1,0,-1,1,0,-1,1,1,1,0,1,0,0,0,1,1,1,0,1,1,0,1,1,0,1,1,1,0]] + for mesh, expect1, expect2 in zip(meshes, vec1_expected, vec2_expected): + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + + print(mesh.entity_orientations) + print(V.cell_node_list) + u = TestFunction(V) + res1= assemble(interpolate(vec1, V)) + for i in range(len(res1.dat.data)): + print(f"{i}: expected: {expect1[i]}, actual: {res1.dat.data[i]}") + assert all(expect1 * res1.dat.data >= 0) + + res2= assemble(interpolate(vec2, V)) + for i in range(len(res2.dat.data)): + print(f"{i}: expected: {expect2[i]}, actual: {res2.dat.data[i]}") + assert all(expect2 * res2.dat.data >= 0) + + +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1)]) +# (construct_nd2_for_fiat, "N1curl", 2), +# (construct_nd2, "N1curl", 2)]) +def test_surface_const_nd(elem_gen, elem_code, deg): cell = polygon(3) elem = construct_nd(cell) ones = as_vector((0,1)) - for n in range(2, 3): + for n in range(1, 6): mesh = UnitSquareMesh(n, n) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, "N1curl", 1) + print(V.cell_node_list) normal = FacetNormal(mesh) ones1 = interpolate(ones, V) res1= assemble(dot(ones1, normal) * ds) @@ -165,21 +179,21 @@ def test_surface_vec(): rt_elem = construct_rt(cell) nd_elem = construct_nd(cell) - for n in range(2, 3): + for n in range(1, 6): mesh = UnitSquareMesh(n, n) x, y = SpatialCoordinate(mesh) normal = FacetNormal(mesh) test_vec = as_vector((-y,x)) - #if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - # V = FunctionSpace(mesh, rt_elem.to_ufl()) - # vec1 = interpolate(test_vec, V) - # res1 = assemble(dot(vec1, normal) * ds) - #else: - # V = FunctionSpace(mesh, "RT", 1) - # vec2 = interpolate(test_vec, V) - # res1 = assemble(dot(vec2, normal) * ds) - #print(f"div {n}: {res1}") + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, rt_elem.to_ufl()) + vec1 = interpolate(test_vec, V) + res1 = assemble(dot(vec1, normal) * ds) + else: + V = FunctionSpace(mesh, "RT", 1) + vec2 = interpolate(test_vec, V) + res1 = assemble(dot(vec2, normal) * ds) + print(f"div {n}: {res1}") if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, nd_elem.to_ufl()) From a6971395dd00564639393ed090e7c9f09c092611 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 28 Oct 2025 17:19:11 +0000 Subject: [PATCH 04/12] work on parameterisation kernel --- fuse/__init__.py | 2 +- fuse/cells.py | 8 ++--- fuse/dof.py | 59 +++++++++++++++++++----------- fuse/triples.py | 1 - fuse/utils.py | 17 +++++---- test/test_orientations.py | 75 ++++++++++++--------------------------- 6 files changed, 76 insertions(+), 86 deletions(-) diff --git a/fuse/__init__.py b/fuse/__init__.py index c0ffec1..dfd083c 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,7 +1,7 @@ from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, tri_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group -from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, PolynomialKernel, ComponentKernel +from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, PolynomialKernel, ComponentKernel, ParameterisationKernel from fuse.triples import ElementTriple, DOFGenerator, immerse from fuse.traces import TrH1, TrGrad, TrHess, TrHCurl, TrHDiv from fuse.tensor_products import tensor_product diff --git a/fuse/cells.py b/fuse/cells.py index b9b03db..16782f4 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -677,7 +677,8 @@ def generate_facet_parameterisation(self, facet_num): res = 0 for d in range(facet_dim + 1): res += coeffs[d] * symbols[d] - return res + return res, symbols[1:] + def plot(self, show=True, plain=False, ax=None, filename=None): """ for now into 2 dimensional space """ @@ -865,10 +866,7 @@ def __call__(self, *x): if hasattr(self.attachment, '__iter__'): res = [] for attach_comp in self.attachment: - 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.append(sympy_to_numpy(attach_comp, syms, x)) return tuple(res) return sympy_to_numpy(self.attachment, syms, x) return x diff --git a/fuse/dof.py b/fuse/dof.py index e51ccdf..3f4f7cb 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -106,6 +106,8 @@ def add_entity(self, entity): def permute(self, g): if g.perm.is_Identity: return self + if self.orientation is not None: + g = self.orientation * g res = L2Pairing() if self.entity: @@ -284,53 +286,62 @@ class ParameterisationKernel(BaseKernel): def __init__(self, g=None): self.g = g self.fn = None + self.syms = None super(ParameterisationKernel, self).__init__() - def add_context(self, cell, entity, fn=None): - new_cls = ParameterisationKernel() - new_cls.cell = cell - new_cls.entity = entity - new_cls.fn = cell.generate_facet_parameterisation(entity.id) + def add_context(self, fn=None, syms=None): + new_cls = ParameterisationKernel(self.g) + new_cls.fn = fn + new_cls.syms = syms return new_cls def __repr__(self): - return f"[{self.comp}]" + if self.g is None or self.g.perm.array_form == [0, 1]: + fn = self.fn + elif self.g.perm.array_form == [1, 0]: + fn = 1 - self.fn + else: + raise NotImplementedError("Cannot parameterise over non-edge subentities.") + return str(fn) def degree(self, interpolant_degree): return interpolant_degree def permute(self, g): + if self.g: + g = self.g * g new_cls = ParameterisationKernel(g) - if self.cell is not None: - new_cls.add_context(self.cell, self.entity) - return new_cls + return new_cls.add_context(self.fn, self.syms) def __call__(self, *args): if self.fn is None: - raise ArgumentError("Function not defined") + raise ValueError("Function not defined") - if self.g is not None: - self.fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) + if self.g.perm.array_form == [1, 0]: + fn = 1 - self.fn + elif self.g.perm.array_form == [0, 1]: + fn = self.fn + else: + raise NotImplementedError("Cannot parameterise over non-edge subentities.") - res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) + res = sympy_to_numpy(fn, self.syms, args[:len(self.syms)]) if not hasattr(res, '__iter__'): - return [res] + res = [res] + #res = np.ones_like(res) return res def tabulate(self, Qpts, attachment=None): - - return np.array([1 for _ in Qpts]).astype(np.float64) return np.array([self(*pt) for pt in Qpts]).astype(np.float64) def _to_dict(self): - o_dict = {"comp": self.comp} + o_dict = {"fn": self.fn} return o_dict def dict_id(self): return "ParameterisationKernel" def _from_dict(obj_dict): - return ParameterisationKernel(obj_dict["comp"]) + return ParameterisationKernel(obj_dict["fn"]) class DOF(): @@ -352,8 +363,6 @@ def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=N self.generation = generation if entity is not None: self.pairing = self.pairing.add_entity(entity) - if self.trace_entity is not None and self.cell is not None: - self.kernel = kernel.add_context(cell, entity) def __call__(self, g): new_generation = self.generation.copy() @@ -379,6 +388,11 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non if self.sub_id is None and generator_id is not None: self.sub_id = generator_id + if self.immersed and isinstance(self.kernel, ParameterisationKernel): + fn, syms = self.cell.generate_facet_parameterisation(self.trace_entity.id) + self.kernel = self.kernel.add_context(fn, syms) + + def convert_to_fiat(self, ref_el, interpolant_degree): return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree) @@ -417,7 +431,7 @@ def eval(self, fn, pullback=True): return self.pairing(self.kernel, attached_fn, self.cell) def tabulate(self, Qpts): - immersion = self.target_space.tabulate(Qpts, self.trace_entity, self.g) + immersion = self.target_space.tabulate(Qpts, self.pairing.entity, self.g) res = self.kernel.tabulate(Qpts, self.attachment) return immersion*res @@ -425,6 +439,9 @@ def __call__(self, g): index_trace = self.cell.d_entities_ids(self.trace_entity.dim()).index(self.trace_entity.id) permuted_e, permuted_g = self.cell.permute_entities(g, self.trace_entity.dim())[index_trace] new_trace_entity = self.cell.get_node(permuted_e).orient(permuted_g) + if self.immersed and isinstance(self.kernel, ParameterisationKernel): + fn, syms = self.cell.generate_facet_parameterisation(new_trace_entity.id) + self.kernel = self.kernel.add_context(fn, syms) new_attach = lambda *x: g(self.attachment(*x)) return ImmersedDOF(self.pairing.permute(permuted_g), self.kernel.permute(permuted_g), new_trace_entity, new_attach, self.target_space, g, self.triple, self.generation, self.sub_id, self.cell) diff --git a/fuse/triples.py b/fuse/triples.py index 309c399..acd6cfc 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -127,7 +127,6 @@ def to_fiat(self): if not pure_perm: self.matrices = mat_perms self.reverse_dof_perms() - #self.orient_mat_perms() form_degree = 1 if self.spaces[0].set_shape else 0 # TODO: Change this when Dense case in Firedrake diff --git a/fuse/utils.py b/fuse/utils.py index 63d1e48..3b3e069 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -17,7 +17,8 @@ def fold_reduce(func_list, *prev): def sympy_to_numpy(array, symbols, values): """ - Convert a sympy array to a numpy array. + TODO: rename this function + Evaluate symbols at values, then convert to numpy if all have been replaced :param: array: sympy array :param: symbols: array of symbols contained in the sympy exprs @@ -27,13 +28,17 @@ def sympy_to_numpy(array, symbols, values): is greater than 1 dimension to remove extra dimensions """ substituted = array.subs({symbols[i]: values[i] for i in range(len(values))}) - nparray = np.array(substituted).astype(np.float64) - if len(nparray.shape) > 1: - return nparray.squeeze() + 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) == 0: - return nparray.item() + if len(nparray.shape) > 1: + return nparray.squeeze() + + if len(nparray.shape) == 0: + return nparray.item() + else: + nparray = substituted return nparray diff --git a/test/test_orientations.py b/test/test_orientations.py index 177aff4..f05edce 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -29,13 +29,10 @@ def construct_nd2(tri=None): x = sp.Symbol("x") y = sp.Symbol("y") - #xs = [DOF(L2Pairing(), PolynomialKernel((x + 1)/2, symbols=[x])), - # DOF(L2Pairing(), PolynomialKernel(1))] - xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), - DOF(L2Pairing(), PolynomialKernel((-np.sqrt(2)/np.sqrt(3))*(6*x + 3), symbols=[x]))] + xs = [DOF(L2Pairing(), ParameterisationKernel())] - dofs = DOFGenerator(xs, S1, S2) + dofs = DOFGenerator(xs, S2, S2) int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) @@ -49,7 +46,6 @@ def construct_nd2(tri=None): Pk = PolynomialSpace(deg - 1) M = sp.Matrix([[y, -x]]) nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M - print(nd.degree()) ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) @@ -97,54 +93,20 @@ def construct_nd2_for_fiat(tri=None): ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) return ned -@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1)]) -# (construct_nd2_for_fiat, "N1curl", 2), -# (construct_nd2, "N1curl", 2)]) -def test_interpolation_values(elem_gen, elem_code,deg): - # this test may not actually make sense proceed with caution - #os.environ["FIREDRAKE_USE_FUSE"] = "1" - cell = polygon(3) - elem = elem_gen(cell) - print() - meshes = [UnitTriangleMesh(0), UnitSquareMesh(1,1), UnitSquareMesh(2,2), UnitSquareMesh(3,3)] - vec1 = as_vector((0, 1)) - vec1_expected = [[0, 1, 1], [0, -1, -1, -1, 0], [1, -1, 0,0,1,0,1,1,0,-1,-1,1,0,-1,0,1], [1,-1,0,0,1,0,1,1,0,1,1,-1,0,0,-1,-1,1,1,-1,-1,0,0,-1,-1,0,1,-1,0,1,0,-1,0,1]] - vec2 = as_vector((1, 0)) - vec2_expected = [[1, -1, 0], [1, 0, 1, 0, -1], [0,1,1,1,0,1,0,-1,1,0,1,0,1,1,1,0], [0,1,1,1,0,1,0,-1,1,0,-1,1,1,1,0,1,0,0,0,1,1,1,0,1,1,0,1,1,0,1,1,1,0]] - for mesh, expect1, expect2 in zip(meshes, vec1_expected, vec2_expected): - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V = FunctionSpace(mesh, elem.to_ufl()) - else: - V = FunctionSpace(mesh, elem_code, deg) - - print(mesh.entity_orientations) - print(V.cell_node_list) - u = TestFunction(V) - res1= assemble(interpolate(vec1, V)) - for i in range(len(res1.dat.data)): - print(f"{i}: expected: {expect1[i]}, actual: {res1.dat.data[i]}") - assert all(expect1 * res1.dat.data >= 0) - - res2= assemble(interpolate(vec2, V)) - for i in range(len(res2.dat.data)): - print(f"{i}: expected: {expect2[i]}, actual: {res2.dat.data[i]}") - assert all(expect2 * res2.dat.data >= 0) - - -@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1)]) -# (construct_nd2_for_fiat, "N1curl", 2), -# (construct_nd2, "N1curl", 2)]) +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), + (construct_nd2_for_fiat, "N1curl", 2), + (construct_nd2, "N1curl", 2)]) def test_surface_const_nd(elem_gen, elem_code, deg): cell = polygon(3) - elem = construct_nd(cell) + elem = elem_gen(cell) ones = as_vector((0,1)) - for n in range(1, 6): + for n in range(1, 2): mesh = UnitSquareMesh(n, n) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: - V = FunctionSpace(mesh, "N1curl", 1) + V = FunctionSpace(mesh, elem_code, deg) print(V.cell_node_list) normal = FacetNormal(mesh) ones1 = interpolate(ones, V) @@ -274,7 +236,6 @@ def test_interpolation(elem_gen, elem_code, deg): def test_create_fiat_nd(): cell = polygon(3) nd = construct_nd2(cell) - breakpoint() nd_fv = construct_nd2_for_fiat(cell) ref_el = cell.to_fiat() sd = ref_el.get_spatial_dimension() @@ -291,12 +252,22 @@ def test_create_fiat_nd(): print(d.convert_to_fiat(ref_el, deg).pt_dict) print("fuse") - for d in nd.generate(): - print(d.convert_to_fiat(ref_el, deg).pt_dict) - #print(d) + dofs = nd.generate() + e = cell.edges()[0] + g = S3.add_cell(cell).members()[2] + print(g) + nodes = [d.convert_to_fiat(ref_el, deg) for d in dofs] + new_nodes = [d(g).convert_to_fiat(ref_el, deg) if d.trace_entity == e else d.convert_to_fiat(ref_el, deg) for d in dofs] + for i in range(len(new_nodes)): + print(f"{dofs[i]}: {nodes[i].pt_dict}") + print(f"{dofs[i]}: {new_nodes[i].pt_dict}") + #for g in S2.add_cell(cell).members(): + # print(d(g)) + # print(f"{g} {d(g).convert_to_fiat(ref_el, deg).pt_dict}") + print() my_elem = nd.to_fiat() - my_elem = nd_fv.to_fiat() - breakpoint() + #my_elem = nd_fv.to_fiat() + #breakpoint() From e76d1d6824b3ba2826a209de08edab19b959752d Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 30 Oct 2025 12:30:00 +0000 Subject: [PATCH 05/12] working on dof permutations --- fuse/cells.py | 2 ++ fuse/dof.py | 7 +++-- fuse/triples.py | 6 +++++ test/test_dofs.py | 8 ++++++ test/test_orientations.py | 55 +++++++++++++++++++-------------------- 5 files changed, 48 insertions(+), 30 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 16782f4..5840c2e 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -803,6 +803,8 @@ def cell_attachment(self, dst): def orient(self, o): """ Orientation node is always labelled with -1 """ + #if self.oriented: + # o = self.oriented * o oriented_point = copy.deepcopy(self) top_level_node = oriented_point.d_entities_ids( oriented_point.dimension)[0] diff --git a/fuse/dof.py b/fuse/dof.py index 3f4f7cb..5f9152c 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -106,12 +106,12 @@ def add_entity(self, entity): def permute(self, g): if g.perm.is_Identity: return self - if self.orientation is not None: - g = self.orientation * g res = L2Pairing() if self.entity: res.entity = self.entity.orient(g) + #if self.orientation is not None: + # g = self.orientation * g res.orientation = g return res @@ -140,6 +140,8 @@ def convert_to_fiat(self, ref_el, dof, interpolant_degree): return functional def __repr__(self): + if self.orientation is not None: + return "integral_{}({})({{kernel}} * {{fn}}) dx) ".format(str(self.orientation), str(self.entity)) return "integral_{}({{kernel}} * {{fn}}) dx) ".format(str(self.entity)) def dict_id(self): @@ -431,6 +433,7 @@ def eval(self, fn, pullback=True): return self.pairing(self.kernel, attached_fn, self.cell) def tabulate(self, Qpts): + # modify this to take reference space q pts immersion = self.target_space.tabulate(Qpts, self.pairing.entity, self.g) res = self.kernel.tabulate(Qpts, self.attachment) return immersion*res diff --git a/fuse/triples.py b/fuse/triples.py index acd6cfc..bab3630 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -123,6 +123,7 @@ def to_fiat(self): nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) counter += 1 self.matrices_by_entity = self.make_entity_dense_matrices(ref_el, entity_ids, nodes, poly_set) + self.matrices_by_entity = self.dummy_function(ref_el, entity_ids, nodes, poly_set) # TODO remove mat_perms, entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) if not pure_perm: self.matrices = mat_perms @@ -136,6 +137,9 @@ def to_fiat(self): dual = DualSet(nodes, ref_el, entity_ids) return CiarletElement(poly_set, dual, degree, form_degree) + def dummy_function(self, ref_el, entity_ids, nodes, poly_set): + return self.matrices_by_entity + def to_tikz(self, show=True, scale=3): """Generates tikz code for the element diagram @@ -273,6 +277,8 @@ def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): new_nodes = [d(g).convert_to_fiat(ref_el, degree) if d.trace_entity == e else d.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)[np.ix_(dof_ids, dof_ids)] + #if dim == 1 and permuted_g.perm.array_form == [1,0]: + # breakpoint() return res_dict def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): diff --git a/test/test_dofs.py b/test/test_dofs.py index ed81804..d007fd7 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -173,3 +173,11 @@ def test_edge_parametrisation(): ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) for n in ned.generate(): print(n) + + from test_orientations import construct_nd2_for_fiat + + ned_fiat = construct_nd2_for_fiat(tri) + + print("fiat") + for n in ned_fiat.generate(): + print(n) diff --git a/test/test_orientations.py b/test/test_orientations.py index f05edce..43b9c36 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -9,17 +9,14 @@ #with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): def dummy_dof_perms(cls, *args, **kwargs): # return -1s of right shape here - oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) - for key1, val1 in oriented_mats_by_entity.items(): - for key2, val2 in oriented_mats_by_entity[key1].items(): - for key3, val3 in oriented_mats_by_entity[key1][key2].items(): - if key1 == 2: - oriented_mats_by_entity[key1][key2][key3] = 1 * np.identity(val3.shape[0]) - #oriented_mats_by_entity[key1][key2][key3][0] = key1*100 + key2*10 + key3 - if key3 == 5: - - oriented_mats_by_entity[key1][key2][key3] = 100 * np.identity(val3.shape[0]) - return oriented_mats_by_entity, False, None + #oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) + oriented_mats_by_entity = cls.make_entity_dense_matrices(*args) + #for key1, val1 in oriented_mats_by_entity.items(): + # for key2, val2 in oriented_mats_by_entity[key1].items(): + # for key3, val3 in oriented_mats_by_entity[key1][key2].items(): + # if key1 == 1 and key3 == 1: + # oriented_mats_by_entity[key1][key2][key3] = np.array([[0, 1],[1,0]]) + return oriented_mats_by_entity def construct_nd2(tri=None): if tri is None: @@ -97,23 +94,25 @@ def construct_nd2_for_fiat(tri=None): (construct_nd2_for_fiat, "N1curl", 2), (construct_nd2, "N1curl", 2)]) def test_surface_const_nd(elem_gen, elem_code, deg): - cell = polygon(3) - elem = elem_gen(cell) - ones = as_vector((0,1)) - - for n in range(1, 2): - mesh = UnitSquareMesh(n, n) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V = FunctionSpace(mesh, elem.to_ufl()) - else: - V = FunctionSpace(mesh, elem_code, deg) - print(V.cell_node_list) - normal = FacetNormal(mesh) - ones1 = interpolate(ones, V) - res1= assemble(dot(ones1, normal) * ds) - - print(f"{n}: {res1}") - assert np.allclose(res1, 0) + + with mock.patch.object(ElementTriple, 'dummy_function', new=dummy_dof_perms): + cell = polygon(3) + elem = elem_gen(cell) + ones = as_vector((0,1)) + + for n in range(2, 3): + mesh = UnitSquareMesh(n, n) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + print(V.cell_node_list) + normal = FacetNormal(mesh) + ones1 = interpolate(ones, V) + res1= assemble(dot(ones1, normal) * ds) + + print(f"{n}: {res1}") + assert np.allclose(res1, 0) def test_surface_const_rt(): From 29e3fcb5eef9e79a203c4ae282efa371367375b8 Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Wed, 5 Nov 2025 10:53:36 +0000 Subject: [PATCH 06/12] Make to fiat follow the definition closer (#40) kernels evaluated on immersed entity cells provide quadrature kernels modify weights and components points are then immersed if needed --- .github/workflows/setup_repos.sh | 18 ++--- .github/workflows/test.yml | 16 +---- Makefile | 12 ++-- conftest.py | 4 +- docs/source/conf.py | 9 ++- fuse/__init__.py | 1 - fuse/cells.py | 24 +++---- fuse/dof.py | 114 +++++++++++++++---------------- fuse/groups.py | 4 +- fuse/spaces/polynomial_spaces.py | 4 +- fuse/traces.py | 10 +-- fuse/triples.py | 75 ++++++++++---------- fuse/utils.py | 2 +- test/test_2d_examples_docs.py | 2 +- test/test_convert_to_fiat.py | 6 +- test/test_dofs.py | 104 +++++++++++++++++----------- test/test_orientations.py | 97 +++++++++++++------------- 17 files changed, 255 insertions(+), 247 deletions(-) diff --git a/.github/workflows/setup_repos.sh b/.github/workflows/setup_repos.sh index f986e1f..ea5d1fe 100644 --- a/.github/workflows/setup_repos.sh +++ b/.github/workflows/setup_repos.sh @@ -16,12 +16,12 @@ git checkout indiamai/integrate_fuse git status python3 -m pip install --break-system-packages -e . -/usr/bin/git config --global --add safe.directory ~ -cd ~ -git clone https://github.com/firedrakeproject/ufl.git -/usr/bin/git config --global --add safe.directory ~/ufl -cd ufl -git fetch -git checkout indiamai/integrate-fuse -git status -python3 -m pip install --break-system-packages -e . \ No newline at end of file +#/usr/bin/git config --global --add safe.directory ~ +#cd ~ +#git clone https://github.com/firedrakeproject/ufl.git +#/usr/bin/git config --global --add safe.directory ~/ufl +#cd ufl +#git fetch +#git checkout indiamai/integrate-fuse +#git status +#python3 -m pip install --break-system-packages -e . diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index a265da9..88b3d0f 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -14,7 +14,8 @@ jobs: # Run on the Github hosted runner runs-on: ubuntu-latest container: - image: firedrakeproject/firedrake-vanilla-default:latest + #image: firedrakeproject/firedrake-vanilla-default:latest + image: firedrakeproject/firedrake-vanilla-default:dev-main # Steps represent a sequence of tasks that will be executed as # part of the jobs steps: @@ -44,17 +45,6 @@ jobs: git checkout indiamai/integrate_fuse git status python3 -m pip install --break-system-packages -e . - - name: Checkout correct UFL branch - run: | - /usr/bin/git config --global --add safe.directory ~ - cd ~ - git clone https://github.com/firedrakeproject/ufl.git - /usr/bin/git config --global --add safe.directory ~/ufl - cd ufl - git fetch - git checkout indiamai/integrate-fuse - git status - python3 -m pip install --break-system-packages -e . - name: Run tests run: | pip list @@ -104,4 +94,4 @@ jobs: message: ${{ env.total }}% minColorRange: 50 maxColorRange: 90 - valColorRange: ${{ env.total }} \ No newline at end of file + valColorRange: ${{ env.total }} diff --git a/Makefile b/Makefile index 969fd02..734a8d3 100644 --- a/Makefile +++ b/Makefile @@ -19,12 +19,12 @@ lint: test_examples: @echo " Running examples" - @python3 -m pytest test/test_2d_examples_docs.py - @python3 -m pytest test/test_3d_examples_docs.py + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_2d_examples_docs.py + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_3d_examples_docs.py tests: @echo " Running all tests" - @python3 -m coverage run -p -m pytest -rx test + @FIREDRAKE_USE_FUSE=1 python3 -m coverage run -p -m pytest -rx test coverage: @python3 -m coverage combine @@ -34,13 +34,13 @@ coverage: test_cells: @echo " Running all cell comparison tests" @firedrake-clean - @python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect0] + @FIREDRAKE_USE_FUSE=1 python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect0] @firedrake-clean - @python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect1] + @FIREDRAKE_USE_FUSE=1 python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect1] clean: @(cd docs/ && make clean) prepush: lint tests @rm .coverage.* - make clean docs \ No newline at end of file + make clean docs diff --git a/conftest.py b/conftest.py index c8c4e00..acf59ca 100644 --- a/conftest.py +++ b/conftest.py @@ -1,9 +1,7 @@ -import pytest - def pytest_addoption(parser): parser.addoption( "--run-cleared", action="store_true", default=False, help="Run tests that require a cleared cache", - ) \ No newline at end of file + ) diff --git a/docs/source/conf.py b/docs/source/conf.py index 154171e..a8781a8 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -87,8 +87,11 @@ # so a file named "default.css" will overwrite the builtin "default.css". html_static_path = ['_static'] # html_css_files = ["additional.css"] + + def setup(app): - app.add_css_file('additional.css') + app.add_css_file('additional.css') + html_theme_options = { 'navbar_links': [ @@ -109,11 +112,11 @@ def setup(app): 'globaltoc_depth': 2, } -html_sidebars = {'api/*': ['localtoc.html'], +html_sidebars = {'api/*': ['localtoc.html'], 'api': ['localtoc.html'], 'about': [], 'install': [], 'index': ['localtoc.html'], 'manual': ['localtoc.html'], 'manual/*': ['localtoc.html'], - '_generated/*': ['custom-sidebar.html'], } \ No newline at end of file + '_generated/*': ['custom-sidebar.html'], } diff --git a/fuse/__init__.py b/fuse/__init__.py index dfd083c..7480548 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,4 +1,3 @@ - from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, tri_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, PolynomialKernel, ComponentKernel, ParameterisationKernel diff --git a/fuse/cells.py b/fuse/cells.py index 5840c2e..9cf7544 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -12,6 +12,7 @@ from sympy.combinatorics.named_groups import SymmetricGroup from fuse.utils import sympy_to_numpy, fold_reduce, numpy_to_str_tuple, orientation_value from FIAT.reference_element import Simplex, TensorProductCell as FiatTensorProductCell, Hypercube +from FIAT.quadrature_schemes import create_quadrature from ufl.cell import Cell, TensorProductCell from functools import cache @@ -656,9 +657,9 @@ def to_tikz(self, show=True, scale=3): return "\n".join(tikz_commands) return tikz_commands - def generate_facet_parameterisation(self, facet_num): - #facet = self.d_entities(self.dimension - 1)[facet_num] + raise NotImplementedError("Facet Parameterisation can be expressed using polynomials") + # facet = self.d_entities(self.dimension - 1)[facet_num] facet = self.get_node(facet_num) facet_dim = facet.dimension if facet_dim != self.dimension - 1: @@ -667,11 +668,9 @@ def generate_facet_parameterisation(self, facet_num): raise NotImplementedError("Facet parameterisation is not implemented for dimensions greater than 1") verts = facet.vertices() v_coords = np.array([self.get_node(v.id, return_coords=True) for v in verts]) - midpoint = np.average(v_coords, axis=0) - stacked = np.c_[np.ones((self.dimension,)), v_coords[:, 0].reshape(self.dimension,1)] + stacked = np.c_[np.ones((self.dimension,)), v_coords[:, 0].reshape(self.dimension, 1)] b = np.array([0, 1]) coeffs = np.linalg.solve(stacked, b) - symbol_names = ["x", "y", "z"] symbols = [1] + [sp.Symbol(symbol_names[d]) for d in range(facet_dim)] res = 0 @@ -679,7 +678,6 @@ def generate_facet_parameterisation(self, facet_num): res += coeffs[d] * symbols[d] return res, symbols[1:] - def plot(self, show=True, plain=False, ax=None, filename=None): """ for now into 2 dimensional space """ if self.dimension == 3: @@ -790,10 +788,12 @@ def attachment(self, source, dst): assert all(np.isclose(val, vals[0]).all() for val in vals) return lambda *x: fold_reduce(attachments[0], *x) - - - - + + def quadrature(self, degree): + fiat_el = self.to_fiat() + Q = create_quadrature(fiat_el, degree) + pts, wts = Q.get_points(), Q.get_weights() + return pts, wts def cell_attachment(self, dst): if not isinstance(dst, int): @@ -803,7 +803,7 @@ def cell_attachment(self, dst): def orient(self, o): """ Orientation node is always labelled with -1 """ - #if self.oriented: + # if self.oriented: # o = self.oriented * o oriented_point = copy.deepcopy(self) top_level_node = oriented_point.d_entities_ids( @@ -1103,7 +1103,7 @@ def __init__(self, cell, name=None): super(CellComplexToUFL, self).__init__(name) def to_fiat(self): - return self.cell_complex.to_fiat(name=self.cellname()) + return self.cell_complex.to_fiat(name=self.cellname) def __repr__(self): return super(CellComplexToUFL, self).__repr__() diff --git a/fuse/dof.py b/fuse/dof.py index 5f9152c..dbb847f 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, IntegralMoment, FrobeniusIntegralMoment +from FIAT.functional import PointEvaluation, IntegralMoment, FrobeniusIntegralMoment, Functional from fuse.utils import sympy_to_numpy import numpy as np import sympy as sp @@ -36,6 +36,7 @@ def __call__(self, kernel, v, cell): return v(*kernel.pt) def convert_to_fiat(self, ref_el, dof, interpolant_deg): + raise NotImplementedError("this should be outdated") pt = dof.eval(FuseFunction(lambda *x: x)) # pt1 = dof.tabulate([[1]]) return PointEvaluation(ref_el, pt) @@ -73,27 +74,11 @@ def __init__(self): super(L2Pairing, self).__init__() def __call__(self, kernel, v, cell): - # TODO get degree of v - # if cell == self.entity: - # ref_el = self.entity.to_fiat() - # # print("evaluating", kernel, v, "on", self.entity) - # Q = create_quadrature(self.entity.to_fiat(), 5) - # # need quadrature here too - therefore need the information from the triple. - # else: - # ref_el = cell.to_fiat() - # ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] - # entity_ref = ref_el.construct_subelement(self.entity.dim()) - # entity = ref_el.construct_subelement(self.entity.dim(), ent_id, self.orientation) - # Q_ref = create_quadrature(entity, 5) - # Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref, self.orientation) - Q = create_quadrature(self.entity.to_fiat(), 5) - - def kernel_dot(x): - return np.dot(kernel(*x), v(*x)) - - return Q.integrate(kernel_dot) + Qpts, Qwts = self.entity.quadrature(5) + return sum([wt*np.dot(kernel(*pt), v(*pt)) for pt, wt in zip(Qpts, Qwts)]) def tabulate(self): + raise NotImplementedError("this should be outdated") pass def add_entity(self, entity): @@ -110,12 +95,13 @@ def permute(self, g): res = L2Pairing() if self.entity: res.entity = self.entity.orient(g) - #if self.orientation is not None: + # if self.orientation is not None: # g = self.orientation * g res.orientation = g return res def convert_to_fiat(self, ref_el, dof, interpolant_degree): + raise NotImplementedError("this should be outdated") total_deg = dof.kernel.degree(interpolant_degree) ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] # entity_ref = ref_el.construct_subelement(self.entity.dim()) @@ -128,8 +114,8 @@ def convert_to_fiat(self, ref_el, dof, interpolant_degree): Jdet = Q.jacobian_determinant() # Jdet = pseudo_determinant(J) qpts, _ = Q.get_points(), Q.get_weights() - # (np.sqrt(2)) * - f_at_qpts = dof.tabulate(qpts).T / Jdet + # (np.sqrt(2)) * + f_at_qpts = dof.tabulate(qpts).T / Jdet comp = None if hasattr(dof.kernel, "comp"): # TODO Value shape should be a variable. @@ -194,10 +180,8 @@ def permute(self, g): def __call__(self, *args): return self.pt - def tabulate(self, Qpts, attachment=None): - #if attachment: - # return np.array([attachment(*self.pt) for _ in Qpts]).astype(np.float64) - return np.array([self.pt for _ in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts): + return np.array([self.pt for _ in Qpts]).astype(np.float64), np.ones_like(Qwts), [[tuple()] for pt in Qpts] def _to_dict(self): o_dict = {"pt": self.pt} @@ -224,13 +208,12 @@ def __repr__(self): def degree(self, interpolant_degree): if len(self.fn.free_symbols) == 0: - return interpolant_degree + return interpolant_degree return self.fn.as_poly().total_degree() * interpolant_degree def permute(self, g): - return self - #new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) - #return PolynomialKernel(new_fn, symbols=self.syms) + new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) + return PolynomialKernel(new_fn, symbols=self.syms) def __call__(self, *args): res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) @@ -238,8 +221,8 @@ def __call__(self, *args): return [res] return res - def tabulate(self, Qpts, attachment=None): - return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts): + return Qpts, np.array([wt*self(*pt) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(len(pt) + 1)] for pt in Qpts] def _to_dict(self): o_dict = {"fn": self.fn} @@ -251,6 +234,7 @@ def dict_id(self): def _from_dict(obj_dict): return PolynomialKernel(obj_dict["fn"]) + class ComponentKernel(BaseKernel): def __init__(self, comp): @@ -258,7 +242,7 @@ def __init__(self, comp): super(ComponentKernel, self).__init__() def __repr__(self): - return f"[{self.comp}]" + return f"[{self.comp}]" def degree(self, interpolant_degree): return interpolant_degree @@ -269,9 +253,9 @@ def permute(self, g): def __call__(self, *args): return args[self.comp] - def tabulate(self, Qpts, attachment=None): - return np.array([1 for _ in Qpts]).astype(np.float64) - #return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts): + return Qpts, Qwts, [[self.comp] for pt in Qpts] + # return Qpts, np.array([self(*pt) for pt in Qpts]).astype(np.float64) def _to_dict(self): o_dict = {"comp": self.comp} @@ -283,9 +267,11 @@ def dict_id(self): def _from_dict(obj_dict): return ComponentKernel(obj_dict["comp"]) + class ParameterisationKernel(BaseKernel): def __init__(self, g=None): + raise NotImplementedError("Parameterisation kernel not needed, use polynomial") self.g = g self.fn = None self.syms = None @@ -293,7 +279,7 @@ def __init__(self, g=None): def add_context(self, fn=None, syms=None): new_cls = ParameterisationKernel(self.g) - new_cls.fn = fn + new_cls.fn = fn new_cls.syms = syms return new_cls @@ -329,11 +315,10 @@ def __call__(self, *args): res = sympy_to_numpy(fn, self.syms, args[:len(self.syms)]) if not hasattr(res, '__iter__'): res = [res] - #res = np.ones_like(res) return res - def tabulate(self, Qpts, attachment=None): - return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts): + return Qpts, np.array([wt*self(*pt) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(len(pt) + 1)] for pt in Qpts] def _to_dict(self): o_dict = {"fn": self.fn} @@ -345,13 +330,14 @@ def dict_id(self): def _from_dict(obj_dict): return ParameterisationKernel(obj_dict["fn"]) + class DOF(): def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=None, g=None, immersed=False, generation=None, sub_id=None, cell=None): self.pairing = pairing self.kernel = kernel self.immersed = immersed - self.trace_entity = entity + self.cell_defined_on = entity self.attachment = attachment self.target_space = target_space self.g = g @@ -368,7 +354,7 @@ def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=N def __call__(self, g): new_generation = self.generation.copy() - return DOF(self.pairing.permute(g), self.kernel.permute(g), self.trace_entity, self.attachment, self.target_space, g, self.immersed, new_generation, self.sub_id, self.cell) + return DOF(self.pairing.permute(g), self.kernel.permute(g), self.cell_defined_on, self.attachment, self.target_space, g, self.immersed, new_generation, self.sub_id, self.cell) def eval(self, fn, pullback=True): return self.pairing(self.kernel, fn, self.cell) @@ -380,8 +366,9 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non # For some of these, we only want to store the first instance of each self.generation[cell.dim()] = dof_gen self.cell = cell - if self.trace_entity is None: + if self.cell_defined_on is None: self.trace_entity = cell + self.cell_defined_on = cell self.pairing = self.pairing.add_entity(cell) if self.target_space is None: self.target_space = space @@ -391,12 +378,25 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non self.sub_id = generator_id if self.immersed and isinstance(self.kernel, ParameterisationKernel): - fn, syms = self.cell.generate_facet_parameterisation(self.trace_entity.id) + fn, syms = self.cell.generate_facet_parameterisation(self.cell_defined_on.id) self.kernel = self.kernel.add_context(fn, syms) - - def convert_to_fiat(self, ref_el, interpolant_degree): - return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree) + 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)) + + def to_quadrature(self, arg_degree): + Qpts, Qwts = self.cell_defined_on.quadrature(arg_degree) + Qwts = Qwts.reshape(Qwts.shape + (1,)) + pts, wts, comps = self.kernel.evaluate(Qpts, Qwts) + if self.immersed: + # need to compute jacobian from attachment. + pts = [self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts] + immersion = self.target_space.tabulate(wts, self.pairing.entity, self.g) + wts = np.outer(wts, immersion) + # pt dict is { pt: (weight, component)} + pt_dict = {tuple(pt): [(w, c) for w, c in zip(wt, cp)] for pt, wt, cp in zip(pts, wts, comps)} + return pt_dict def __repr__(self, fn="v"): return str(self.pairing).format(fn=fn, kernel=self.kernel) @@ -428,29 +428,29 @@ def eval(self, fn, pullback=True): attached_fn = fn.attach(self.attachment) if pullback: - attached_fn = self.target_space(attached_fn, self.trace_entity, self.g) + attached_fn = self.target_space(attached_fn, self.cell_defined_on, self.g) return self.pairing(self.kernel, attached_fn, self.cell) def tabulate(self, Qpts): # modify this to take reference space q pts immersion = self.target_space.tabulate(Qpts, self.pairing.entity, self.g) - res = self.kernel.tabulate(Qpts, self.attachment) + res, _ = self.kernel.tabulate(Qpts, self.attachment) return immersion*res def __call__(self, g): - index_trace = self.cell.d_entities_ids(self.trace_entity.dim()).index(self.trace_entity.id) - permuted_e, permuted_g = self.cell.permute_entities(g, self.trace_entity.dim())[index_trace] - new_trace_entity = self.cell.get_node(permuted_e).orient(permuted_g) + index_trace = self.cell.d_entities_ids(self.cell_defined_on.dim()).index(self.cell_defined_on.id) + permuted_e, permuted_g = self.cell.permute_entities(g, self.cell_defined_on.dim())[index_trace] + new_cell_defined_on = self.cell.get_node(permuted_e).orient(permuted_g) if self.immersed and isinstance(self.kernel, ParameterisationKernel): - fn, syms = self.cell.generate_facet_parameterisation(new_trace_entity.id) + fn, syms = self.cell.generate_facet_parameterisation(new_cell_defined_on.id) self.kernel = self.kernel.add_context(fn, syms) new_attach = lambda *x: g(self.attachment(*x)) - return ImmersedDOF(self.pairing.permute(permuted_g), self.kernel.permute(permuted_g), new_trace_entity, + return ImmersedDOF(self.pairing.permute(permuted_g), self.kernel.permute(permuted_g), new_cell_defined_on, new_attach, self.target_space, g, self.triple, self.generation, self.sub_id, self.cell) def __repr__(self): - fn = "tr_{1}_{0}(v)".format(str(self.trace_entity), str(self.target_space)) + fn = "tr_{1}_{0}(v)".format(str(self.cell_defined_on), str(self.target_space)) return super(ImmersedDOF, self).__repr__(fn) def immerse(self, entity, attachment, trace, g): @@ -505,5 +505,3 @@ def _to_dict(self): def dict_id(self): return "Function" - - diff --git a/fuse/groups.py b/fuse/groups.py index 14a5a41..3ddc77c 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -165,12 +165,11 @@ def get_member(self, perm): return m raise ValueError("Permutation not a member of group") - def get_member_by_val(self, val): for m in self.members(): if m.numeric_rep() == val: return m - + raise ValueError("Value does not represent a group member") def compute_num_reps(self, base_val=0): @@ -285,7 +284,6 @@ def get_member(self, perm): if m.perm == perm: return m raise ValueError("Permutation not a member of group") - # def compute_reps(self, g, path, remaining_members): # # breadth first search to find generator representations of all members diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 56bbce1..422cf05 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -147,8 +147,8 @@ 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_mat(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/traces.py b/fuse/traces.py index 2f77606..e283dfb 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -15,7 +15,7 @@ def __call__(self, trace_entity, g): def plot(self, ax, coord, trace_entity, g, **kwargs): raise NotImplementedError("Trace uninstanitated") - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity, g): raise NotImplementedError("Tabulation uninstantiated") def _to_dict(self): @@ -54,8 +54,8 @@ def plot(self, ax, coord, trace_entity, g, **kwargs): def to_tikz(self, coord, trace_entity, g, scale, color="black"): return f"\\filldraw[{color}] {numpy_to_str_tuple(coord, scale)} circle (2pt) node[anchor = south] {{}};" - def tabulate(self, Qpts, trace_entity, g): - return np.ones_like(Qpts) + def tabulate(self, Qwts, trace_entity, g): + return Qwts def __repr__(self): return "H1" @@ -80,7 +80,7 @@ def plot(self, ax, coord, trace_entity, g, **kwargs): vec = self.tabulate([], trace_entity, g).squeeze() ax.quiver(*coord, *vec, **kwargs) - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity, g): entityBasis = np.array(trace_entity.basis_vectors()) cellEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) basis = np.matmul(entityBasis, cellEntityBasis) @@ -117,7 +117,7 @@ def apply(*x): return tuple(result) return apply - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity, g): tangent = np.array(trace_entity.basis_vectors()) subEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) result = np.matmul(tangent, subEntityBasis) diff --git a/fuse/triples.py b/fuse/triples.py index bab3630..fb4ab21 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -77,16 +77,16 @@ def degree(self): def get_dof_info(self, dof, tikz=True): colours = {False: {0: "b", 1: "r", 2: "g", 3: "b"}, True: {0: "blue", 1: "red", 2: "green", 3: "black"}} - if dof.trace_entity.dimension == 0: - center = self.cell.cell_attachment(dof.trace_entity.id)() - elif dof.trace_entity.dimension == 1: - center = self.cell.cell_attachment(dof.trace_entity.id)(0) - elif dof.trace_entity.dimension == 2: - center = self.cell.cell_attachment(dof.trace_entity.id)(0, 0) + if dof.cell_defined_on.dimension == 0: + center = self.cell.cell_attachment(dof.cell_defined_on.id)() + elif dof.cell_defined_on.dimension == 1: + center = self.cell.cell_attachment(dof.cell_defined_on.id)(0) + elif dof.cell_defined_on.dimension == 2: + center = self.cell.cell_attachment(dof.cell_defined_on.id)(0, 0) else: center = list(sum(np.array(self.cell.vertices(return_coords=True)))) - return center, colours[tikz][dof.trace_entity.dimension] + return center, colours[tikz][dof.cell_defined_on.dimension] def get_value_shape(self): # TODO Shape should be specificed somewhere else probably @@ -102,6 +102,7 @@ def to_fiat(self): ref_el = self.cell.to_fiat() dofs = self.generate() degree = self.spaces[0].degree() + value_shape = self.get_value_shape() entity_ids = {} entity_perms = {} nodes = [] @@ -118,16 +119,18 @@ def to_fiat(self): for entity in entities: 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) - nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) + if entity[1] == dofs[i].cell_defined_on.id - min_ids[dim]: + entity_ids[dim][dofs[i].cell_defined_on.id - min_ids[dim]].append(counter) + nodes.append(dofs[i].convert_to_fiat(ref_el, degree, value_shape)) counter += 1 self.matrices_by_entity = self.make_entity_dense_matrices(ref_el, entity_ids, nodes, poly_set) - self.matrices_by_entity = self.dummy_function(ref_el, entity_ids, nodes, poly_set) # TODO remove + self.matrices_by_entity = self.dummy_function(ref_el, entity_ids, nodes, poly_set) # TODO remove mat_perms, entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) if not pure_perm: self.matrices = mat_perms self.reverse_dof_perms() + else: + self.matrices = entity_perms form_degree = 1 if self.spaces[0].set_shape else 0 # TODO: Change this when Dense case in Firedrake @@ -157,12 +160,12 @@ def to_tikz(self, show=True, scale=3): if isinstance(dof.pairing, DeltaPairing): coord = dof.eval(identity, pullback=False) if isinstance(dof.target_space, Trace): - tikz_commands += [dof.target_space.to_tikz(coord, dof.trace_entity, dof.g, scale, color)] + tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, dof.g, scale, color)] else: tikz_commands += [f"\\filldraw[{color}] {numpy_to_str_tuple(coord, scale)} circle (2pt) node[anchor = south] {{}};"] elif isinstance(dof.pairing, L2Pairing): coord = center - tikz_commands += [dof.target_space.to_tikz(coord, dof.trace_entity, dof.g, scale, color)] + tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, dof.g, scale, color)] if show: tikz_commands += ['\\end{tikzpicture}'] return "\n".join(tikz_commands) @@ -190,7 +193,7 @@ def plot(self, filename="temp.png"): if len(coord) == 1: coord = (coord[0], 0) if isinstance(dof.target_space, Trace): - dof.target_space.plot(ax, coord, dof.trace_entity, dof.g, color=color) + dof.target_space.plot(ax, coord, dof.cell_defined_on, dof.g, color=color) else: ax.scatter(*coord, color=color) ax.text(*coord, dof.id) @@ -212,12 +215,12 @@ def plot(self, filename="temp.png"): if isinstance(dof.pairing, DeltaPairing): coord = dof.eval(identity, pullback=False) if isinstance(dof.target_space, Trace): - dof.target_space.plot(ax, coord, dof.trace_entity, dof.g, color=color) + dof.target_space.plot(ax, coord, dof.cell_defined_on, dof.g, color=color) else: ax.scatter(*coord, color=color) elif isinstance(dof.pairing, L2Pairing): coord = center - dof.target_space.plot(ax, center, dof.trace_entity, dof.g, color=color, length=0.2) + dof.target_space.plot(ax, center, dof.cell_defined_on, dof.g, color=color, length=0.2) ax.text(*coord, dof.id) plt.axis('off') ax.get_xaxis().set_visible(False) @@ -251,7 +254,7 @@ def compute_dense_matrix(self, ref_el, entity_ids, nodes, poly_set): def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): degree = self.spaces[0].degree() min_ids = self.cell.get_starter_ids() - nodes = [d.convert_to_fiat(ref_el, degree) for d in self.generate()] + nodes = [d.convert_to_fiat(ref_el, degree, self.get_value_shape()) for d in self.generate()] sub_ents = [] res_dict = {} for d in range(0, self.cell.dim() + 1): @@ -261,7 +264,7 @@ def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): dim = e.dim() e_id = e.id - min_ids[dim] res_dict[dim][e_id] = {} - dof_ids = [d.id for d in self.generate() if d.trace_entity == e] + dof_ids = [d.id for d in self.generate() if d.cell_defined_on == e] res_dict[dim][e_id][0] = np.eye(len(dof_ids)) original_V, original_basis = self.compute_dense_matrix(ref_el, entity_ids, nodes, poly_set) @@ -274,11 +277,11 @@ def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): elif g.perm.is_Identity: res_dict[dim][e_id][val] = np.eye(len(dof_ids)) else: - new_nodes = [d(g).convert_to_fiat(ref_el, degree) if d.trace_entity == e else d.convert_to_fiat(ref_el, degree) for d in self.generate()] + new_nodes = [d(g).convert_to_fiat(ref_el, degree, self.get_value_shape()) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, degree, self.get_value_shape()) 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)[np.ix_(dof_ids, dof_ids)] - #if dim == 1 and permuted_g.perm.array_form == [1,0]: - # breakpoint() + # if dim == 1 and permuted_g.perm.array_form == [1,0]: + # breakpoint() return res_dict def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): @@ -294,7 +297,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 = [d(g).convert_to_fiat(ref_el, degree, self.get_value_shape()) 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 @@ -310,8 +313,8 @@ 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.trace_entity.dim() - sub_dict = entity_associations[sub_dim][d.trace_entity.id - min_ids[sub_dim]] + 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]): dof_gen = str(d.generation[dim]) @@ -356,9 +359,9 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): dofs = self.generate() min_ids = self.cell.get_starter_ids() entity_associations, pure_perm, sub_pure_perm = self._entity_associations(dofs) - #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), None, pure_perm + # 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), None, pure_perm # return self.matrices_by_entity, None, pure_perm oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(dofs) @@ -378,8 +381,9 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): # (dof_gen, ent_dofs) total_ent_dof_ids += [ed.id for ed in ent_dofs if ed.id not in total_ent_dof_ids] dof_gen_class = ent_dofs[0].generation - if not len(dof_gen_class[dim].g2.members()) == 1: + if not len(dof_gen_class[dim].g2.members()) == 1 and dim == min(dof_gen_class.keys()): # if DOFs on entity are not perms, get the matrix + # only get this if they are defined on the current dimension oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = self.matrices_by_entity[dim][e_id][val] elif g.perm.is_Identity or (pure_perm and len(ent_dofs_ids) == 1): oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = np.eye(len(ent_dofs_ids)) @@ -393,11 +397,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() else: # TODO what if an orientation is not in G1 - #warnings.warn("FUSE: should probably be using equivalent members of g in g1 for this") - #sub_mat = g.matrix_form() - #oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() - - #raise NotImplementedError(f"Orientation {g} is not in group {dof_gen_class[dim].g1.members()}") + warnings.warn("FUSE: orientation case not covered") + # sub_mat = g.matrix_form() + # oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() + # raise NotImplementedError(f"Orientation {g} is not in group {dof_gen_class[dim].g1.members()}") pass if len(dof_gen_class.keys()) == 2 and dim == self.cell.dim(): @@ -433,6 +436,7 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): return oriented_mats_by_entity, None, False def orient_mat_perms(self): + raise NotImplementedError("This should not be necessary") min_ids = self.cell.get_starter_ids() entity_orientations = compare_topologies(ufc_cell(self.cell.to_ufl().cellname()).get_topology(), self.cell.get_topology()) num_ents = 0 @@ -440,8 +444,7 @@ def orient_mat_perms(self): ents = self.cell.d_entities(dim) for e in ents: e_id = e.id - min_ids[dim] - members = e.group.members()#and dim < self.cell.dim() - if entity_orientations[num_ents + e_id] != 0 : + if entity_orientations[num_ents + e_id] != 0: modifier = self.matrices[dim][e_id][entity_orientations[num_ents+e_id]] reverse_modifier_val = (~e.group.get_member_by_val(entity_orientations[num_ents+e_id])).numeric_rep() reverse_modifier = self.matrices[dim][e_id][reverse_modifier_val] @@ -535,7 +538,7 @@ def make_entity_ids(self): entity_ids[dim] = {i: [] for i in top[dim]} for i in range(len(dofs)): - entity = dofs[i].trace_entity + entity = dofs[i].cell_defined_on dim = entity.dim() entity_ids[dim][entity.id - min_ids[dim]].append(i) return entity_ids diff --git a/fuse/utils.py b/fuse/utils.py index 3b3e069..fcfa3d9 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -34,7 +34,7 @@ def sympy_to_numpy(array, symbols, values): if len(nparray.shape) > 1: return nparray.squeeze() - + if len(nparray.shape) == 0: return nparray.item() else: diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index de6be11..8b64f6e 100644 --- a/test/test_2d_examples_docs.py +++ b/test/test_2d_examples_docs.py @@ -165,7 +165,7 @@ def construct_nd(tri=None): y = sp.Symbol("y") # xs = [DOF(L2Pairing(), PointKernel(edge.basis_vectors()[0]))] - #xs = [DOF(L2Pairing(), PointKernel((np.sqrt(2),)))] + # xs = [DOF(L2Pairing(), PointKernel((np.sqrt(2),)))] xs = [DOF(L2Pairing(), PolynomialKernel((1,)))] dofs = DOFGenerator(xs, S1, S2) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 128a53e..f405009 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -470,7 +470,7 @@ def test_helmholtz_3d(elem_gen, elem_code, deg, conv_rate): def helmholtz_solve(V, mesh): # Define variational problem - dim = mesh.ufl_cell().topological_dimension() + dim = mesh.ufl_cell().topological_dimension u = TrialFunction(V) v = TestFunction(V) f = Function(V) @@ -594,10 +594,10 @@ def test_project_vec(elem_gen, elem_code, deg): mesh = UnitTriangleMesh() U = FunctionSpace(mesh, elem_code, deg) - assert np.allclose(project(U, mesh, as_vector((1,1))), 0, rtol=1e-5) + assert np.allclose(project(U, mesh, as_vector((1, 1))), 0, rtol=1e-5) U = FunctionSpace(mesh, elem.to_ufl()) - assert np.allclose(project(U, mesh, as_vector((1,1))), 0, rtol=1e-5) + assert np.allclose(project(U, mesh, as_vector((1, 1))), 0, rtol=1e-5) @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_dg1_tet, "DG", 1)]) diff --git a/test/test_dofs.py b/test/test_dofs.py index d007fd7..7afbbe6 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -1,5 +1,7 @@ from fuse import * from test_convert_to_fiat import create_cg1, create_dg1, construct_cg3, construct_rt, construct_nd +from test_orientations import construct_nd2 + import sympy as sp import numpy as np @@ -139,45 +141,65 @@ def test_permute_nodes(): for g in cg1.cell.group.members(): print(g, [d(g).convert_to_fiat(cell.to_fiat(), degree).pt_dict for d in cg1.generate()]) -def test_edge_parametrisation(): - tri = polygon(3) - for i in tri.d_entities(1): - print(tri.generate_facet_parameterisation(i.id)) - from fuse.dof import ParameterisationKernel - - - deg = 2 - edge = tri.edges()[0] - x = sp.Symbol("x") - y = sp.Symbol("y") - - xs = [DOF(L2Pairing(), ParameterisationKernel())] - - - dofs = DOFGenerator(xs, S2, S2) - int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - - xs = [DOF(L2Pairing(), ComponentKernel((0,))), - DOF(L2Pairing(), ComponentKernel((1,)))] - center_dofs = DOFGenerator(xs, S1, S3) - xs = [immerse(tri, int_ned1, TrHCurl)] - tri_dofs = DOFGenerator(xs, C3, S1) - - vec_Pk = PolynomialSpace(deg - 1, set_shape=True) - Pk = PolynomialSpace(deg - 1) - M = sp.Matrix([[y, -x]]) - nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M - - - ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) - for n in ned.generate(): - print(n) - - from test_orientations import construct_nd2_for_fiat - - ned_fiat = construct_nd2_for_fiat(tri) - - print("fiat") - for n in ned_fiat.generate(): - print(n) +# def test_edge_parametrisation(): +# tri = polygon(3) +# for i in tri.d_entities(1): +# print(tri.generate_facet_parameterisation(i.id)) +# from fuse.dof import ParameterisationKernel +# +# deg = 2 +# edge = tri.edges()[0] +# x = sp.Symbol("x") +# y = sp.Symbol("y") +# +# xs = [DOF(L2Pairing(), ParameterisationKernel())] +# +# dofs = DOFGenerator(xs, S2, S2) +# int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) +# +# xs = [DOF(L2Pairing(), ComponentKernel((0,))), +# DOF(L2Pairing(), ComponentKernel((1,)))] +# center_dofs = DOFGenerator(xs, S1, S3) +# xs = [immerse(tri, int_ned1, TrHCurl)] +# tri_dofs = DOFGenerator(xs, C3, S1) +# +# vec_Pk = PolynomialSpace(deg - 1, set_shape=True) +# Pk = PolynomialSpace(deg - 1) +# M = sp.Matrix([[y, -x]]) +# nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M +# +# ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) +# for n in ned.generate(): +# print(n) +# +# from test_orientations import construct_nd2_for_fiat +# +# ned_fiat = construct_nd2_for_fiat(tri) +# +# print("fiat") +# for n in ned_fiat.generate(): +# print(n) + + +def test_generate_quadrature(): + cell = polygon(3) + degree = 2 + # elem = create_dg1(cell) + # elem = create_cg1(cell) + # elem = construct_nd(cell) + elem = construct_nd2(cell) + # elem = construct_nd2_for_fiat(cell) + from FIAT.nedelec import Nedelec + fiat_elem = Nedelec(cell.to_fiat(), degree) + # from FIAT.lagrange import Lagrange + # fiat_elem = Lagrange(cell.to_fiat(), degree) + degree = elem.spaces[0].degree() + print(degree) + for d in fiat_elem.dual_basis(): + print("fiat", d.pt_dict) + print() + for d in elem.generate(): + print("fuse", d.to_quadrature(degree)) + + elem.to_fiat() diff --git a/test/test_orientations.py b/test/test_orientations.py index 43b9c36..77a7cb3 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -3,21 +3,23 @@ from firedrake import * from fuse import * import sympy as sp -from test_convert_to_fiat import create_cg1, helmholtz_solve, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg0, create_dg1 +from test_convert_to_fiat import create_cg1, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg1 import os -#with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): + +# with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): def dummy_dof_perms(cls, *args, **kwargs): # return -1s of right shape here - #oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) + # oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) oriented_mats_by_entity = cls.make_entity_dense_matrices(*args) - #for key1, val1 in oriented_mats_by_entity.items(): + # for key1, val1 in oriented_mats_by_entity.items(): # for key2, val2 in oriented_mats_by_entity[key1].items(): # for key3, val3 in oriented_mats_by_entity[key1][key2].items(): # if key1 == 1 and key3 == 1: - # oriented_mats_by_entity[key1][key2][key3] = np.array([[0, 1],[1,0]]) + # oriented_mats_by_entity[key1][key2][key3] = np.array([[0, 1],[1,0]]) return oriented_mats_by_entity + def construct_nd2(tri=None): if tri is None: tri = polygon(3) @@ -26,28 +28,26 @@ def construct_nd2(tri=None): x = sp.Symbol("x") y = sp.Symbol("y") - xs = [DOF(L2Pairing(), ParameterisationKernel())] - + xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] dofs = DOFGenerator(xs, S2, S2) int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - xs = [DOF(L2Pairing(), ComponentKernel((0,))), DOF(L2Pairing(), ComponentKernel((1,)))] center_dofs = DOFGenerator(xs, S1, S3) - xs = [immerse(tri, int_ned1, TrHCurl)] + xs = [immerse(tri, int_ned1, TrHCurl)] tri_dofs = DOFGenerator(xs, C3, S1) vec_Pk = PolynomialSpace(deg - 1, set_shape=True) Pk = PolynomialSpace(deg - 1) M = sp.Matrix([[y, -x]]) - nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M - + nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) return ned + def construct_nd2_for_fiat(tri=None): if tri is None: tri = polygon(3) @@ -59,7 +59,6 @@ def construct_nd2_for_fiat(tri=None): xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), DOF(L2Pairing(), PolynomialKernel((3*x*np.sqrt(2)/np.sqrt(3)), symbols=[x]))] - dofs = DOFGenerator(xs, S1, S2) int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) @@ -76,7 +75,7 @@ def construct_nd2_for_fiat(tri=None): xs = [DOF(L2Pairing(), ComponentKernel((0,))), DOF(L2Pairing(), ComponentKernel((1,)))] center_dofs = DOFGenerator(xs, S1, S3) - xs = [immerse(tri, int_ned2, TrHCurl, node=0)] + xs = [immerse(tri, int_ned2, TrHCurl, node=0)] xs += [immerse(tri, int_ned1, TrHCurl, node=1)] xs += [immerse(tri, int_ned3, TrHCurl, node=2)] tri_dofs = DOFGenerator(xs, S1, S1) @@ -84,23 +83,21 @@ def construct_nd2_for_fiat(tri=None): vec_Pk = PolynomialSpace(deg - 1, set_shape=True) Pk = PolynomialSpace(deg - 1) M = sp.Matrix([[y, -x]]) - nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M - + nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs]) return ned + @pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), - (construct_nd2_for_fiat, "N1curl", 2), (construct_nd2, "N1curl", 2)]) def test_surface_const_nd(elem_gen, elem_code, deg): - with mock.patch.object(ElementTriple, 'dummy_function', new=dummy_dof_perms): cell = polygon(3) elem = elem_gen(cell) - ones = as_vector((0,1)) + ones = as_vector((0, 1)) - for n in range(2, 3): + for n in range(2, 6): mesh = UnitSquareMesh(n, n) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) @@ -109,8 +106,8 @@ def test_surface_const_nd(elem_gen, elem_code, deg): print(V.cell_node_list) normal = FacetNormal(mesh) ones1 = interpolate(ones, V) - res1= assemble(dot(ones1, normal) * ds) - + res1 = assemble(dot(ones1, normal) * ds) + print(f"{n}: {res1}") assert np.allclose(res1, 0) @@ -118,7 +115,7 @@ def test_surface_const_nd(elem_gen, elem_code, deg): def test_surface_const_rt(): cell = polygon(3) elem = construct_rt(cell) - ones = as_vector((1,0)) + ones = as_vector((1, 0)) for n in range(1, 6): mesh = UnitSquareMesh(n, n) @@ -129,12 +126,13 @@ def test_surface_const_rt(): V = FunctionSpace(mesh, "RT", 1) normal = FacetNormal(mesh) ones1 = interpolate(ones, V) - res1= assemble(dot(ones1, normal) * ds) - + res1 = assemble(dot(ones1, normal) * ds) + print(f"{n}: {res1}") assert np.allclose(res1, 0) +@pytest.mark.xfail(reason="orientations") def test_surface_vec(): cell = polygon(3) rt_elem = construct_rt(cell) @@ -145,7 +143,7 @@ def test_surface_vec(): mesh = UnitSquareMesh(n, n) x, y = SpatialCoordinate(mesh) normal = FacetNormal(mesh) - test_vec = as_vector((-y,x)) + test_vec = as_vector((-y, x)) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, rt_elem.to_ufl()) vec1 = interpolate(test_vec, V) @@ -155,13 +153,13 @@ def test_surface_vec(): vec2 = interpolate(test_vec, V) res1 = assemble(dot(vec2, normal) * ds) print(f"div {n}: {res1}") - + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, nd_elem.to_ufl()) print(V.cell_node_list) vec1 = interpolate(test_vec, V) - #V2 = FunctionSpace(mesh, create_dg0(cell).to_ufl()) - #u = TestFunction(V2) + # V2 = FunctionSpace(mesh, create_dg0(cell).to_ufl()) + # u = TestFunction(V2) res2 = assemble(dot(vec1, normal) * ds) else: V = FunctionSpace(mesh, "N1curl", 1) @@ -170,13 +168,13 @@ def test_surface_vec(): res2 = assemble(dot(vec1, normal) * ds) print(f"curl {n}: {res2}") - #assert np.allclose(0, res1) + assert np.allclose(0, res1) assert np.allclose(0, res2) -def test_interpolate_vs_project(V): +def interpolate_vs_project(V): mesh = V.mesh() - dim = mesh.geometric_dimension() + dim = mesh.geometric_dimension if dim == 2: x, y = SpatialCoordinate(mesh) elif dim == 3: @@ -207,14 +205,15 @@ def test_interpolate_vs_project(V): def test_degree2_interpolation(): cell = polygon(3) - elem = construct_nd2_for_fiat(cell) - mesh = UnitSquareMesh(1,1) + elem = construct_nd2(cell) + mesh = UnitSquareMesh(1, 1) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, "N1curl", 2) - test_interpolate_vs_project(V) + interpolate_vs_project(V) + @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), (create_cg1, "CG", 1), @@ -224,20 +223,20 @@ def test_degree2_interpolation(): def test_interpolation(elem_gen, elem_code, deg): cell = polygon(3) elem = elem_gen(cell) - mesh = UnitSquareMesh(1,1) + mesh = UnitSquareMesh(1, 1) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, elem_code, deg) - test_interpolate_vs_project(V) + interpolate_vs_project(V) + def test_create_fiat_nd(): cell = polygon(3) nd = construct_nd2(cell) nd_fv = construct_nd2_for_fiat(cell) ref_el = cell.to_fiat() - sd = ref_el.get_spatial_dimension() deg = 2 from FIAT.nedelec import Nedelec @@ -246,7 +245,7 @@ def test_create_fiat_nd(): for f in fiat_elem.dual_basis(): print(f.pt_dict) - print("fiat: fuse's version") + print("fiat: fuse's version") for d in nd_fv.generate(): print(d.convert_to_fiat(ref_el, deg).pt_dict) @@ -255,18 +254,16 @@ def test_create_fiat_nd(): e = cell.edges()[0] g = S3.add_cell(cell).members()[2] print(g) - nodes = [d.convert_to_fiat(ref_el, deg) for d in dofs] - new_nodes = [d(g).convert_to_fiat(ref_el, deg) if d.trace_entity == e else d.convert_to_fiat(ref_el, deg) for d in dofs] + nodes = [d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] + new_nodes = [d(g).convert_to_fiat(ref_el, deg, (2,)) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] for i in range(len(new_nodes)): - print(f"{dofs[i]}: {nodes[i].pt_dict}") - print(f"{dofs[i]}: {new_nodes[i].pt_dict}") - #for g in S2.add_cell(cell).members(): - # print(d(g)) - # print(f"{g} {d(g).convert_to_fiat(ref_el, deg).pt_dict}") + print(f"{nodes[i].pt_dict}") + # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") + # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") + # for g in S2.add_cell(cell).members(): + # print(d(g)) + # print(f"{g} {d(g).convert_to_fiat(ref_el, deg).pt_dict}") print() - my_elem = nd.to_fiat() - #my_elem = nd_fv.to_fiat() - #breakpoint() - - + nd.to_fiat() + nd_fv.to_fiat() From 48a6612c782b0be8630f1d6d8d914b69ad7b3aa2 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 13 Nov 2025 18:53:58 +0000 Subject: [PATCH 07/12] work on orientation of dofs - restructuring, still needs a big clean up --- fuse/__init__.py | 2 +- fuse/cells.py | 11 ++- fuse/dof.py | 171 ++++++++++++-------------------------- fuse/groups.py | 12 +-- fuse/traces.py | 55 ++++++------ fuse/triples.py | 25 +++--- test/test_dofs.py | 23 ++++- test/test_orientations.py | 88 +++++++++++++------- 8 files changed, 189 insertions(+), 198 deletions(-) diff --git a/fuse/__init__.py b/fuse/__init__.py index 7480548..80c5970 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,6 +1,6 @@ from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, tri_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group -from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, PolynomialKernel, ComponentKernel, ParameterisationKernel +from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, PolynomialKernel, ComponentKernel from fuse.triples import ElementTriple, DOFGenerator, immerse from fuse.traces import TrH1, TrGrad, TrHess, TrHCurl, TrHDiv from fuse.tensor_products import tensor_product diff --git a/fuse/cells.py b/fuse/cells.py index 9cf7544..036b35e 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -624,7 +624,11 @@ def basis_vectors(self, return_coords=True, entity=None): if self.dimension == 0: # return [[] raise ValueError("Dimension 0 entities cannot have Basis Vectors") - top_level_node = self_levels[0][0] + if self.oriented: + # ordered_vertices() handles the orientation so we want to drop the orientation node + top_level_node = self_levels[1][0] + else: + top_level_node = self_levels[0][0] v_0 = vertices[0] if return_coords: v_0_coords = self.attachment(top_level_node, v_0)() @@ -803,8 +807,8 @@ def cell_attachment(self, dst): def orient(self, o): """ Orientation node is always labelled with -1 """ - # if self.oriented: - # o = self.oriented * o + if self.oriented: + o = self.oriented * o oriented_point = copy.deepcopy(self) top_level_node = oriented_point.d_entities_ids( oriented_point.dimension)[0] @@ -912,7 +916,6 @@ def get_spatial_dimension(self): def get_sub_entities(self): self.A.get_sub_entities() self.B.get_sub_entities() - breakpoint() def dimension(self): return tuple(self.A.dimension, self.B.dimension) diff --git a/fuse/dof.py b/fuse/dof.py index dbb847f..556efec 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -35,12 +35,8 @@ def __call__(self, kernel, v, cell): assert isinstance(kernel, PointKernel) return v(*kernel.pt) - def convert_to_fiat(self, ref_el, dof, interpolant_deg): - raise NotImplementedError("this should be outdated") - pt = dof.eval(FuseFunction(lambda *x: x)) - # pt1 = dof.tabulate([[1]]) - return PointEvaluation(ref_el, pt) - # return PointEvaluation(ref_el, tuple(pt1[0])) + def tabulate(self): + return 1 def add_entity(self, entity): res = DeltaPairing() @@ -52,7 +48,8 @@ def add_entity(self, entity): def permute(self, g): res = DeltaPairing() if self.entity: - res.entity = self.entity.orient(g) + res.entity = self.entity + # res.entity = self.entity.orient(g) res.orientation = g return res @@ -72,14 +69,18 @@ class L2Pairing(Pairing): def __init__(self): super(L2Pairing, self).__init__() - def __call__(self, kernel, v, cell): Qpts, Qwts = self.entity.quadrature(5) return sum([wt*np.dot(kernel(*pt), v(*pt)) for pt, wt in zip(Qpts, Qwts)]) def tabulate(self): - raise NotImplementedError("this should be outdated") - pass + + bvs = np.array(self.entity.basis_vectors()) + if self.orientation: + new_bvs = np.array(self.entity.orient(self.orientation).basis_vectors()) + basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) + return basis_change + return np.eye(bvs.shape[0]) def add_entity(self, entity): res = L2Pairing() @@ -94,37 +95,15 @@ def permute(self, g): res = L2Pairing() if self.entity: - res.entity = self.entity.orient(g) - # if self.orientation is not None: - # g = self.orientation * g + res.entity = self.entity + #res.entity = self.entity.orient(g) + #if self.entity.oriented: + # breakpoint() + if self.orientation is not None: + g = self.orientation * g res.orientation = g return res - def convert_to_fiat(self, ref_el, dof, interpolant_degree): - raise NotImplementedError("this should be outdated") - total_deg = dof.kernel.degree(interpolant_degree) - ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] - # entity_ref = ref_el.construct_subelement(self.entity.dim()) - entity = ref_el.construct_subelement(self.entity.dim(), ent_id, self.orientation) - Q_ref = create_quadrature(entity, total_deg) - # pts_ref, wts_ref = Q_ref.get_points(), Q_ref.get_weights() - - # pts, wts, J = map_quadrature(pts_ref, wts_ref, Q_ref.ref_el, entity_ref, jacobian=True) - Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref, self.orientation) - Jdet = Q.jacobian_determinant() - # Jdet = pseudo_determinant(J) - qpts, _ = Q.get_points(), Q.get_weights() - # (np.sqrt(2)) * - f_at_qpts = dof.tabulate(qpts).T / Jdet - comp = None - if hasattr(dof.kernel, "comp"): - # TODO Value shape should be a variable. - comp = dof.kernel.comp - functional = IntegralMoment(ref_el, Q, f_at_qpts, comp=comp, shp=(2,)) - else: - functional = FrobeniusIntegralMoment(ref_el, Q, f_at_qpts) - return functional - def __repr__(self): if self.orientation is not None: return "integral_{}({})({{kernel}} * {{fn}}) dx) ".format(str(self.orientation), str(self.entity)) @@ -180,7 +159,7 @@ def permute(self, g): def __call__(self, *args): return self.pt - def evaluate(self, Qpts, Qwts): + def evaluate(self, Qpts, Qwts, basis_change): return np.array([self.pt for _ in Qpts]).astype(np.float64), np.ones_like(Qwts), [[tuple()] for pt in Qpts] def _to_dict(self): @@ -196,10 +175,11 @@ def _from_dict(obj_dict): class PolynomialKernel(BaseKernel): - def __init__(self, fn, symbols=[]): + def __init__(self, fn, g=None, symbols=[]): if len(symbols) != 0 and not sp.sympify(fn).as_poly(): raise ValueError("Function argument must be able to be interpreted as a sympy polynomial") self.fn = sp.sympify(fn) + self.g = g self.syms = symbols super(PolynomialKernel, self).__init__() @@ -213,16 +193,19 @@ 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))}) - return PolynomialKernel(new_fn, symbols=self.syms) + return PolynomialKernel(new_fn, g=g, symbols=self.syms) def __call__(self, *args): res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) if not hasattr(res, '__iter__'): return [res] + #if self.g: + # print(self.g, self.g(res)) + # return self.g(res) return res - def evaluate(self, Qpts, Qwts): - return Qpts, np.array([wt*self(*pt) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(len(pt) + 1)] for pt in Qpts] + def evaluate(self, Qpts, Qwts, basis_change): + 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(len(pt) + 1)] for pt in Qpts] def _to_dict(self): o_dict = {"fn": self.fn} @@ -238,6 +221,7 @@ def _from_dict(obj_dict): class ComponentKernel(BaseKernel): def __init__(self, comp): + self.comp = comp super(ComponentKernel, self).__init__() @@ -251,9 +235,10 @@ def permute(self, g): return self def __call__(self, *args): - return args[self.comp] + return tuple(args[i] if i in self.comp else 0 for i in range(len(args))) +# return tuple(args[c] for c in self.comp) - def evaluate(self, Qpts, Qwts): + def evaluate(self, Qpts, Qwts, basis_change): return Qpts, Qwts, [[self.comp] for pt in Qpts] # return Qpts, np.array([self(*pt) for pt in Qpts]).astype(np.float64) @@ -268,68 +253,6 @@ def _from_dict(obj_dict): return ComponentKernel(obj_dict["comp"]) -class ParameterisationKernel(BaseKernel): - - def __init__(self, g=None): - raise NotImplementedError("Parameterisation kernel not needed, use polynomial") - self.g = g - self.fn = None - self.syms = None - super(ParameterisationKernel, self).__init__() - - def add_context(self, fn=None, syms=None): - new_cls = ParameterisationKernel(self.g) - new_cls.fn = fn - new_cls.syms = syms - return new_cls - - def __repr__(self): - if self.g is None or self.g.perm.array_form == [0, 1]: - fn = self.fn - elif self.g.perm.array_form == [1, 0]: - fn = 1 - self.fn - else: - raise NotImplementedError("Cannot parameterise over non-edge subentities.") - return str(fn) - - def degree(self, interpolant_degree): - return interpolant_degree - - def permute(self, g): - if self.g: - g = self.g * g - new_cls = ParameterisationKernel(g) - return new_cls.add_context(self.fn, self.syms) - - def __call__(self, *args): - if self.fn is None: - raise ValueError("Function not defined") - - if self.g.perm.array_form == [1, 0]: - fn = 1 - self.fn - elif self.g.perm.array_form == [0, 1]: - fn = self.fn - else: - raise NotImplementedError("Cannot parameterise over non-edge subentities.") - - res = sympy_to_numpy(fn, self.syms, args[:len(self.syms)]) - if not hasattr(res, '__iter__'): - res = [res] - return res - - def evaluate(self, Qpts, Qwts): - return Qpts, np.array([wt*self(*pt) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(len(pt) + 1)] for pt in Qpts] - - def _to_dict(self): - o_dict = {"fn": self.fn} - return o_dict - - def dict_id(self): - return "ParameterisationKernel" - - def _from_dict(obj_dict): - return ParameterisationKernel(obj_dict["fn"]) - class DOF(): @@ -377,9 +300,9 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non if self.sub_id is None and generator_id is not None: self.sub_id = generator_id - if self.immersed and isinstance(self.kernel, ParameterisationKernel): - fn, syms = self.cell.generate_facet_parameterisation(self.cell_defined_on.id) - self.kernel = self.kernel.add_context(fn, syms) + #if self.immersed and isinstance(self.kernel, ParameterisationKernel): + # fn, syms = self.cell.generate_facet_parameterisation(self.cell_defined_on.id) + # self.kernel = self.kernel.add_context(fn, syms) def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): # TODO deriv dict needs implementing (currently {}) @@ -388,14 +311,25 @@ def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): def to_quadrature(self, arg_degree): Qpts, Qwts = self.cell_defined_on.quadrature(arg_degree) Qwts = Qwts.reshape(Qwts.shape + (1,)) - pts, wts, comps = self.kernel.evaluate(Qpts, Qwts) + bvs = np.array(self.cell_defined_on.basis_vectors()) + + if self.pairing.orientation: + new_bvs = np.array(self.cell_defined_on.orient(self.pairing.orientation).basis_vectors()) + basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) + else: + basis_change = np.eye(bvs.shape[0]) + pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, np.eye(bvs.shape[0])) if self.immersed: # need to compute jacobian from attachment. pts = [self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts] - immersion = self.target_space.tabulate(wts, self.pairing.entity, self.g) + immersion = self.target_space.tabulate(wts, self.pairing.entity) + immersion = np.matmul(basis_change, immersion) wts = np.outer(wts, immersion) + # pt dict is { pt: (weight, component)} pt_dict = {tuple(pt): [(w, c) for w, c in zip(wt, cp)] for pt, wt, cp in zip(pts, wts, comps)} + #if isinstance(self.kernel, ComponentKernel): + # breakpoint() return pt_dict def __repr__(self, fn="v"): @@ -428,23 +362,24 @@ def eval(self, fn, pullback=True): attached_fn = fn.attach(self.attachment) if pullback: - attached_fn = self.target_space(attached_fn, self.cell_defined_on, self.g) + attached_fn = self.target_space(attached_fn, self.cell_defined_on) return self.pairing(self.kernel, attached_fn, self.cell) def tabulate(self, Qpts): # modify this to take reference space q pts - immersion = self.target_space.tabulate(Qpts, self.pairing.entity, self.g) + immersion = self.target_space.tabulate(Qpts, self.pairing.entity) res, _ = self.kernel.tabulate(Qpts, self.attachment) return immersion*res def __call__(self, g): index_trace = self.cell.d_entities_ids(self.cell_defined_on.dim()).index(self.cell_defined_on.id) permuted_e, permuted_g = self.cell.permute_entities(g, self.cell_defined_on.dim())[index_trace] - new_cell_defined_on = self.cell.get_node(permuted_e).orient(permuted_g) - if self.immersed and isinstance(self.kernel, ParameterisationKernel): - fn, syms = self.cell.generate_facet_parameterisation(new_cell_defined_on.id) - self.kernel = self.kernel.add_context(fn, syms) + new_cell_defined_on = self.cell.get_node(permuted_e) + #.orient(permuted_g) + #if self.immersed and isinstance(self.kernel, ParameterisationKernel): + # fn, syms = self.cell.generate_facet_parameterisation(new_cell_defined_on.id) + # self.kernel = self.kernel.add_context(fn, syms) new_attach = lambda *x: g(self.attachment(*x)) return ImmersedDOF(self.pairing.permute(permuted_g), self.kernel.permute(permuted_g), new_cell_defined_on, new_attach, self.target_space, g, self.triple, self.generation, self.sub_id, self.cell) diff --git a/fuse/groups.py b/fuse/groups.py index 3ddc77c..d9a2a5c 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -279,11 +279,13 @@ def transform_between_perms(self, perm1, perm2): assert perm2 in member_perms return ~self.get_member(Permutation(perm1)) * self.get_member(Permutation(perm2)) - def get_member(self, perm): - for m in self.members(): - if m.perm == perm: - return m - raise ValueError("Permutation not a member of group") + #def get_member(self, perm): + # if not isinstance(perm, Permutation): + # perm = Permutation.from_sequence(perm) + # for m in self.members(): + # if m.perm == perm: + # return m + # raise ValueError("Permutation not a member of group") # def compute_reps(self, g, path, remaining_members): # # breadth first search to find generator representations of all members diff --git a/fuse/traces.py b/fuse/traces.py index e283dfb..ff2ca5e 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -9,13 +9,13 @@ class Trace(): def __init__(self, cell): self.domain = cell - def __call__(self, trace_entity, g): + def __call__(self, trace_entity): raise NotImplementedError("Trace uninstanitated") - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): raise NotImplementedError("Trace uninstanitated") - def tabulate(self, Qwts, trace_entity, g): + def tabulate(self, Qwts, trace_entity): raise NotImplementedError("Tabulation uninstantiated") def _to_dict(self): @@ -45,16 +45,16 @@ class TrH1(Trace): def __init__(self, cell): super(TrH1, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): return v - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): ax.scatter(*coord, **kwargs) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): + def to_tikz(self, coord, trace_entity, scale, color="black"): return f"\\filldraw[{color}] {numpy_to_str_tuple(coord, scale)} circle (2pt) node[anchor = south] {{}};" - def tabulate(self, Qwts, trace_entity, g): + def tabulate(self, Qwts, trace_entity): return Qwts def __repr__(self): @@ -66,21 +66,21 @@ class TrHDiv(Trace): def __init__(self, cell): super(TrHDiv, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): def apply(*x): - result = np.dot(self.tabulate(None, trace_entity, g), np.array(v(*x)).squeeze()) + result = np.dot(self.tabulate(None, trace_entity), np.array(v(*x)).squeeze()) if isinstance(result, np.float64): # todo: might always be a float return (result,) return tuple(result) return apply - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): # plot dofs of the type associated with this space - vec = self.tabulate([], trace_entity, g).squeeze() + vec = self.tabulate([], trace_entity).squeeze() ax.quiver(*coord, *vec, **kwargs) - def tabulate(self, Qwts, trace_entity, g): + def tabulate(self, Qwts, trace_entity): entityBasis = np.array(trace_entity.basis_vectors()) cellEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) basis = np.matmul(entityBasis, cellEntityBasis) @@ -94,8 +94,8 @@ def tabulate(self, Qwts, trace_entity, g): return result - def to_tikz(self, coord, trace_entity, g, scale, color="black"): - vec = self.tabulate([], trace_entity, g).squeeze() + def to_tikz(self, coord, trace_entity, scale, color="black"): + vec = self.tabulate([], trace_entity).squeeze() end_point = [coord[i] + 0.25*vec[i] for i in range(len(coord))] arw = "-{Stealth[length=3mm, width=2mm]}" return f"\\draw[thick, {color}, {arw}] {numpy_to_str_tuple(coord, scale)} -- {numpy_to_str_tuple(end_point, scale)};" @@ -109,26 +109,26 @@ class TrHCurl(Trace): def __init__(self, cell): super(TrHCurl, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): def apply(*x): - result = np.dot(self.tabulate(None, trace_entity, g), np.array(v(*x)).squeeze()) + result = np.dot(self.tabulate(None, trace_entity), np.array(v(*x)).squeeze()) if isinstance(result, np.float64): return (result,) return tuple(result) return apply - def tabulate(self, Qwts, trace_entity, g): + def tabulate(self, Qwts, trace_entity): tangent = np.array(trace_entity.basis_vectors()) subEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) result = np.matmul(tangent, subEntityBasis) return result - def plot(self, ax, coord, trace_entity, g, **kwargs): - vec = self.tabulate([], trace_entity, g).squeeze() + def plot(self, ax, coord, trace_entity, **kwargs): + vec = self.tabulate([], trace_entity).squeeze() ax.quiver(*coord, *vec, **kwargs) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): - vec = self.tabulate([], trace_entity, g).squeeze() + def to_tikz(self, coord, trace_entity, scale, color="black"): + vec = self.tabulate([], trace_entity).squeeze() end_point = [coord[i] + 0.25*vec[i] for i in range(len(coord))] arw = "-{Stealth[length=3mm, width=2mm]}" return f"\\draw[thick, {color}, {arw}] {numpy_to_str_tuple(coord, scale)} -- {numpy_to_str_tuple(end_point, scale)};" @@ -142,7 +142,7 @@ class TrGrad(Trace): def __init__(self, cell): super(TrGrad, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): # Compute grad v and then dot with tangent rotated according to the group member tangent = np.array(g(np.array(self.domain.basis_vectors())[0])) @@ -159,11 +159,11 @@ def apply(*x): return tuple(result) return apply - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): circle1 = plt.Circle(coord, 0.075, fill=False, **kwargs) ax.add_patch(circle1) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): + def to_tikz(self, coord, trace_entity, scale, color="black"): return f"\\draw[{color}] {numpy_to_str_tuple(coord, scale)} circle (4pt) node[anchor = south] {{}};" def __repr__(self): @@ -175,7 +175,8 @@ class TrHess(Trace): def __init__(self, cell): super(TrHess, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): + raise NotImplementedError("Hessian trace needs reviewing") b0, b1 = self.domain.basis_vectors() tangent0 = np.array(g(b0)) tangent1 = np.array(g(b1)) @@ -192,11 +193,11 @@ def apply(*x): return tuple(result) return apply - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): circle1 = plt.Circle(coord, 0.15, fill=False, **kwargs) ax.add_patch(circle1) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): + def to_tikz(self, coord, trace_entity, scale, color="black"): return f"\\draw[{color}] {numpy_to_str_tuple(coord, scale)} circle (6pt) node[anchor = south] {{}};" def __repr__(self): diff --git a/fuse/triples.py b/fuse/triples.py index fb4ab21..45d2b72 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -126,11 +126,12 @@ def to_fiat(self): self.matrices_by_entity = self.make_entity_dense_matrices(ref_el, entity_ids, nodes, poly_set) self.matrices_by_entity = self.dummy_function(ref_el, entity_ids, nodes, poly_set) # TODO remove mat_perms, entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) - if not pure_perm: - self.matrices = mat_perms - self.reverse_dof_perms() - else: - self.matrices = entity_perms + self.matrices = mat_perms + # if not pure_perm: + # self.matrices = mat_perms + # else: + # self.matrices = entity_perms + self.reverse_dof_perms() form_degree = 1 if self.spaces[0].set_shape else 0 # TODO: Change this when Dense case in Firedrake @@ -160,12 +161,12 @@ def to_tikz(self, show=True, scale=3): if isinstance(dof.pairing, DeltaPairing): coord = dof.eval(identity, pullback=False) if isinstance(dof.target_space, Trace): - tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, dof.g, scale, color)] + tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, scale, color)] else: tikz_commands += [f"\\filldraw[{color}] {numpy_to_str_tuple(coord, scale)} circle (2pt) node[anchor = south] {{}};"] elif isinstance(dof.pairing, L2Pairing): coord = center - tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, dof.g, scale, color)] + tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, scale, color)] if show: tikz_commands += ['\\end{tikzpicture}'] return "\n".join(tikz_commands) @@ -193,7 +194,7 @@ def plot(self, filename="temp.png"): if len(coord) == 1: coord = (coord[0], 0) if isinstance(dof.target_space, Trace): - dof.target_space.plot(ax, coord, dof.cell_defined_on, dof.g, color=color) + dof.target_space.plot(ax, coord, dof.cell_defined_on, color=color) else: ax.scatter(*coord, color=color) ax.text(*coord, dof.id) @@ -215,12 +216,12 @@ def plot(self, filename="temp.png"): if isinstance(dof.pairing, DeltaPairing): coord = dof.eval(identity, pullback=False) if isinstance(dof.target_space, Trace): - dof.target_space.plot(ax, coord, dof.cell_defined_on, dof.g, color=color) + dof.target_space.plot(ax, coord, dof.cell_defined_on, color=color) else: ax.scatter(*coord, color=color) elif isinstance(dof.pairing, L2Pairing): coord = center - dof.target_space.plot(ax, center, dof.cell_defined_on, dof.g, color=color, length=0.2) + dof.target_space.plot(ax, center, dof.cell_defined_on, color=color, length=0.2) ax.text(*coord, dof.id) plt.axis('off') ax.get_xaxis().set_visible(False) @@ -280,8 +281,8 @@ def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): new_nodes = [d(g).convert_to_fiat(ref_el, degree, self.get_value_shape()) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, degree, self.get_value_shape()) 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)[np.ix_(dof_ids, dof_ids)] - # if dim == 1 and permuted_g.perm.array_form == [1,0]: - # breakpoint() + #if dim == 1 and permuted_g.perm.array_form == [1,0]: + # breakpoint() return res_dict def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): diff --git a/test/test_dofs.py b/test/test_dofs.py index 7afbbe6..72c0838 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -77,11 +77,32 @@ def test_permute_nd(): # print(dof) for g in nd.cell.group.members(): - print("g:", g.numeric_rep()) + print("g:",g, g.numeric_rep()) for dof in nd.generate(): print(dof(g).convert_to_fiat(cell.to_fiat(), 0).pt_dict) print(dof, "->", dof(g), "eval, ", dof(g).eval(func)) +def test_permute_nd2(): + cell = polygon(3) + + nd = construct_nd2(cell) + x = sp.Symbol("x") + y = sp.Symbol("y") + func = FuseFunction(sp.Matrix([x, -1/3 + 2*y]), symbols=(x, y)) + + # for dof in nd.generate(): + # print(dof) + + for g in nd.cell.group.members(): + print("g:", g.numeric_rep()) + if g.numeric_rep() == 2: + for i,dof in enumerate(nd.generate()): + if 0 < i < 2: + print(dof(g).convert_to_fiat(cell.to_fiat(), 2, (2,)).pt_dict) + print(dof, "->", dof(g), "eval, ", dof(g).eval(func)) + breakpoint() + + def test_permute_nd_old(): cell = polygon(3) diff --git a/test/test_orientations.py b/test/test_orientations.py index 77a7cb3..448987c 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -33,8 +33,8 @@ def construct_nd2(tri=None): dofs = DOFGenerator(xs, S2, S2) int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - xs = [DOF(L2Pairing(), ComponentKernel((0,))), - DOF(L2Pairing(), ComponentKernel((1,)))] + xs = [DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[0])), + DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[1]))] center_dofs = DOFGenerator(xs, S1, S3) xs = [immerse(tri, int_ned1, TrHCurl)] tri_dofs = DOFGenerator(xs, C3, S1) @@ -92,24 +92,23 @@ def construct_nd2_for_fiat(tri=None): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), (construct_nd2, "N1curl", 2)]) def test_surface_const_nd(elem_gen, elem_code, deg): - with mock.patch.object(ElementTriple, 'dummy_function', new=dummy_dof_perms): - cell = polygon(3) - elem = elem_gen(cell) - ones = as_vector((0, 1)) - - for n in range(2, 6): - mesh = UnitSquareMesh(n, n) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V = FunctionSpace(mesh, elem.to_ufl()) - else: - V = FunctionSpace(mesh, elem_code, deg) - print(V.cell_node_list) - normal = FacetNormal(mesh) - ones1 = interpolate(ones, V) - res1 = assemble(dot(ones1, normal) * ds) + cell = polygon(3) + elem = elem_gen(cell) + ones = as_vector((0, 1)) - print(f"{n}: {res1}") - assert np.allclose(res1, 0) + for n in range(2, 6): + mesh = UnitSquareMesh(n, n) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + print(V.cell_node_list) + normal = FacetNormal(mesh) + ones1 = interpolate(ones, V) + res1 = assemble(dot(ones1, normal) * ds) + + print(f"{n}: {res1}") + assert np.allclose(res1, 0) def test_surface_const_rt(): @@ -132,7 +131,7 @@ def test_surface_const_rt(): assert np.allclose(res1, 0) -@pytest.mark.xfail(reason="orientations") +#@pytest.mark.xfail(reason="orientations") def test_surface_vec(): cell = polygon(3) rt_elem = construct_rt(cell) @@ -245,25 +244,54 @@ def test_create_fiat_nd(): for f in fiat_elem.dual_basis(): print(f.pt_dict) - print("fiat: fuse's version") - for d in nd_fv.generate(): - print(d.convert_to_fiat(ref_el, deg).pt_dict) + #print("fiat: fuse's version") + #for d in nd_fv.generate(): + # print(d.convert_to_fiat(ref_el, deg).pt_dict) print("fuse") dofs = nd.generate() e = cell.edges()[0] - g = S3.add_cell(cell).members()[2] + s3 = S3.add_cell(cell) + g = s3.get_member([0,2,1]) + for d in nd.generate(): + print(d.convert_to_fiat(ref_el, deg, (2,)).pt_dict) print(g) - nodes = [d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] - new_nodes = [d(g).convert_to_fiat(ref_el, deg, (2,)) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] - for i in range(len(new_nodes)): - print(f"{nodes[i].pt_dict}") + for d in nd.generate(): + print(d(g).convert_to_fiat(ref_el, deg, (2,)).pt_dict) + #nodes = [d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] + #new_nodes = [d(g).convert_to_fiat(ref_el, deg, (2,)) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] + #for i in range(len(new_nodes)): + # print(f"{nodes[i].pt_dict}") # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") # for g in S2.add_cell(cell).members(): # print(d(g)) # print(f"{g} {d(g).convert_to_fiat(ref_el, deg).pt_dict}") - print() + nd1= construct_nd(cell) + print("nd1") + for d in nd1.generate(): + print(d.convert_to_fiat(ref_el, deg).pt_dict) + print(g) + for d in nd1.generate(): + print(d(g).convert_to_fiat(ref_el, deg).pt_dict) + nd1.to_fiat() + nd.to_fiat() - nd_fv.to_fiat() + breakpoint() + #nd_fv.to_fiat() + + +def test_int_nd(): + tri = polygon(3) + deg = 2 + edge = tri.edges()[0] + x = sp.Symbol("x") + y = sp.Symbol("y") + + xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] + + dofs = DOFGenerator(xs, S2, S2) + int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) + dofs = int_ned1.generate() + breakpoint() From 0ba61916628c7412178685c108e443c28819d08f Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 14 Nov 2025 18:50:32 +0000 Subject: [PATCH 08/12] reassess transformation of kernel, generating basis pair through reflection --- fuse/cells.py | 2 +- fuse/dof.py | 16 +++-- fuse/groups.py | 8 ++- fuse/triples.py | 7 +- test/test_orientations.py | 145 +++++++++++++++++++++++++++----------- 5 files changed, 126 insertions(+), 52 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 036b35e..17e6fd4 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -95,7 +95,7 @@ def compute_scaled_verts(d, n): :param: n: number of vertices """ if d == 2: - source = np.array([0, 1]) + source = np.array([-np.sqrt(3)/2, -1/2]) rot_coords = [source for i in range(0, n)] rot_mat = np.array([[np.cos((2*np.pi)/n), -np.sin((2*np.pi)/n)], [np.sin((2*np.pi)/n), np.cos((2*np.pi)/n)]]) diff --git a/fuse/dof.py b/fuse/dof.py index 556efec..610304d 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -205,7 +205,9 @@ def __call__(self, *args): return res def evaluate(self, Qpts, Qwts, basis_change): - 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(len(pt) + 1)] for pt in Qpts] + # convert basis change to right format + return Qpts, np.array([wt*np.matmul(self(*pt), basis_change) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(len(pt) + 1)] 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), [[(i,) for i in range(len(pt) + 1)] for pt in Qpts] def _to_dict(self): o_dict = {"fn": self.fn} @@ -311,19 +313,21 @@ def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): def to_quadrature(self, arg_degree): Qpts, Qwts = self.cell_defined_on.quadrature(arg_degree) Qwts = Qwts.reshape(Qwts.shape + (1,)) - bvs = np.array(self.cell_defined_on.basis_vectors()) - if self.pairing.orientation: + if self.cell_defined_on.get_spatial_dimension() > 0 and self.pairing.orientation: + bvs = np.array(self.cell_defined_on.basis_vectors()) new_bvs = np.array(self.cell_defined_on.orient(self.pairing.orientation).basis_vectors()) basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) + print(basis_change) else: - basis_change = np.eye(bvs.shape[0]) - pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, np.eye(bvs.shape[0])) + basis_change = np.eye(self.cell_defined_on.get_spatial_dimension()) + pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change) if self.immersed: # need to compute jacobian from attachment. pts = [self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts] immersion = self.target_space.tabulate(wts, self.pairing.entity) - immersion = np.matmul(basis_change, immersion) + #if basis_change: + # immersion = np.matmul(immersion, basis_change) wts = np.outer(wts, immersion) # pt dict is { pt: (weight, component)} diff --git a/fuse/groups.py b/fuse/groups.py index d9a2a5c..d6e7457 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -78,7 +78,11 @@ def __hash__(self): def __mul__(self, x): assert isinstance(x, GroupMemberRep) - return self.group.get_member(self.perm * x.perm) + max_size = max(self.perm.size, x.perm.size) + larger_group = self.group if len(self.group.members()) >= len(x.group.members()) else x.group + x_perm = Permutation(x.perm, size=max_size) + self_perm = Permutation(self.perm, size=max_size) + return larger_group.get_member(self_perm * x_perm) def __invert__(self): return self.group.get_member(~self.perm) @@ -217,6 +221,8 @@ def __init__(self, base_group, cell=None): if cell is not None: self.cell = cell vertices = cell.vertices(return_coords=True) + verts = cell.ordered_vertices() + vertices = [cell.get_node(v, return_coords=True) for v in verts] self._members = [] counter = 0 diff --git a/fuse/triples.py b/fuse/triples.py index 45d2b72..617e99c 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -122,6 +122,7 @@ def to_fiat(self): if entity[1] == dofs[i].cell_defined_on.id - min_ids[dim]: entity_ids[dim][dofs[i].cell_defined_on.id - min_ids[dim]].append(counter) nodes.append(dofs[i].convert_to_fiat(ref_el, degree, value_shape)) + print(nodes[-1].pt_dict) counter += 1 self.matrices_by_entity = self.make_entity_dense_matrices(ref_el, entity_ids, nodes, poly_set) self.matrices_by_entity = self.dummy_function(ref_el, entity_ids, nodes, poly_set) # TODO remove @@ -429,10 +430,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): # remove immersed DOFs from flat form oriented_mats_overall = oriented_mats_by_entity[dim][0] - for val, mat in oriented_mats_overall.items(): - cell_dofs = entity_ids[dim][0] - flat_by_entity[dim][e_id][val] = perm_matrix_to_perm_array(mat[np.ix_(cell_dofs, cell_dofs)]) if pure_perm and sub_pure_perm: + for val, mat in oriented_mats_overall.items(): + cell_dofs = entity_ids[dim][0] + flat_by_entity[dim][e_id][val] = perm_matrix_to_perm_array(mat[np.ix_(cell_dofs, cell_dofs)]) return oriented_mats_by_entity, flat_by_entity, True return oriented_mats_by_entity, None, False diff --git a/test/test_orientations.py b/test/test_orientations.py index 448987c..566252d 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -7,19 +7,6 @@ import os -# with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): -def dummy_dof_perms(cls, *args, **kwargs): - # return -1s of right shape here - # oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) - oriented_mats_by_entity = cls.make_entity_dense_matrices(*args) - # for key1, val1 in oriented_mats_by_entity.items(): - # for key2, val2 in oriented_mats_by_entity[key1].items(): - # for key3, val3 in oriented_mats_by_entity[key1][key2].items(): - # if key1 == 1 and key3 == 1: - # oriented_mats_by_entity[key1][key2][key3] = np.array([[0, 1],[1,0]]) - return oriented_mats_by_entity - - def construct_nd2(tri=None): if tri is None: tri = polygon(3) @@ -32,10 +19,12 @@ def construct_nd2(tri=None): dofs = DOFGenerator(xs, S2, S2) int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - - xs = [DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[0])), - DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[1]))] - center_dofs = DOFGenerator(xs, S1, S3) + v_2 = np.array(tri.get_node(tri.ordered_vertices()[2], return_coords = True)) + v_1 = np.array(tri.get_node(tri.ordered_vertices()[1], return_coords = True)) + xs = [DOF(L2Pairing(), PolynomialKernel(v_2 - v_1))] + #xs = [DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[0]))] +# DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[1]))] + center_dofs = DOFGenerator(xs, S2, S3) xs = [immerse(tri, int_ned1, TrHCurl)] tri_dofs = DOFGenerator(xs, C3, S1) @@ -131,11 +120,12 @@ def test_surface_const_rt(): assert np.allclose(res1, 0) -#@pytest.mark.xfail(reason="orientations") -def test_surface_vec(): +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), + (construct_nd2, "N1curl", 2)]) +def test_surface_vec(elem_gen,elem_code,deg): cell = polygon(3) rt_elem = construct_rt(cell) - nd_elem = construct_nd(cell) + nd_elem = elem_gen(cell) for n in range(1, 6): @@ -161,7 +151,7 @@ def test_surface_vec(): # u = TestFunction(V2) res2 = assemble(dot(vec1, normal) * ds) else: - V = FunctionSpace(mesh, "N1curl", 1) + V = FunctionSpace(mesh, elem_code, deg) print(V.cell_node_list) vec1 = interpolate(test_vec, V) res2 = assemble(dot(vec1, normal) * ds) @@ -182,36 +172,57 @@ def interpolate_vs_project(V): shape = V.value_shape if dim == 2: if len(shape) == 0: - expression = x + y + exact = Function(FunctionSpace(mesh, 'CG', 5)) + #expression = x + y + expression = x elif len(shape) == 1: - expression = as_vector([x, y]) - elif len(shape) == 2: - expression = as_tensor(([x, y], [x, y])) + exact = Function(VectorFunctionSpace(mesh, 'CG', 5)) + expr = x + y + #expr = cos(x*pi*2)*sin(y*pi*2) + expression = as_vector([expr, expr]) elif dim == 3: if len(shape) == 0: + exact = Function(FunctionSpace(mesh, 'CG', 5)) expression = x + y + z elif len(shape) == 1: + exact = Function(FunctionSpace(mesh, 'CG', 5)) expression = as_vector([x, y, z]) - elif len(shape) == 2: - expression = as_tensor(([x, y, z], [x, y, z], [x, y, z])) - f = assemble(interpolate(expression, V)) + breakpoint() expect = project(expression, V) - print(f.dat.data) - print(expect.dat.data) - assert np.allclose(f.dat.data, expect.dat.data, atol=1e-06) + exact.interpolate(expression) + return sqrt(assemble(inner((expect - exact), (expect - exact)) * dx)), sqrt(assemble(inner((f - exact), (f - exact)) * dx)), -def test_degree2_interpolation(): - cell = polygon(3) - elem = construct_nd2(cell) - mesh = UnitSquareMesh(1, 1) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V = FunctionSpace(mesh, elem.to_ufl()) - else: - V = FunctionSpace(mesh, "N1curl", 2) - interpolate_vs_project(V) + +@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [(construct_nd2, "N1curl", 2, 1.8)]) +def test_convergence(elem_gen,elem_code,deg,conv_rate): + cell = polygon(3) + elem = elem_gen(cell) + scale_range = range(2,6) + diff_proj = [0 for i in scale_range] + diff_inte = [0 for i in scale_range] + for n in scale_range: + mesh = UnitSquareMesh(2**n, 2**n) + + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + diff_proj[n-1], diff_inte[n-1] = interpolate_vs_project(V) + + breakpoint() + print("projection l2 error norms:", diff_proj) + diff_proj = np.array(diff_proj) + conv1 = np.log2(diff_proj[:-1] / diff_proj[1:]) + print("convergence order:", conv1) + assert all([c > conv_rate for c in conv1]) + print("interpolation l2 error norms:", diff_inte) + diff_inte = np.array(diff_inte) + conv1 = np.log2(diff_inte[:-1] / diff_inte[1:]) + print("convergence order:", conv1) + assert all([c > conv_rate for c in conv1]) @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), @@ -228,8 +239,50 @@ def test_interpolation(elem_gen, elem_code, deg): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, elem_code, deg) - interpolate_vs_project(V) + err1, err2 = interpolate_vs_project(V) + assert np.allclose(err1, 0) + assert np.allclose(err2, 0) + +def test_projection_convergence(elem_gen, elem_code, deg, conv_rate): + cell = make_tetrahedron() + elem = elem_gen(cell) + function = lambda x: cos((3/4)*pi*x[0]) + scale_range = range(1, 6) + diff = [0 for i in scale_range] + for i in scale_range: + mesh = UnitSquareMesh(2 ** i, 2 ** i) + x, y = SpatialCoordinate(mesh) + expr = cos(x*pi*2)*sin(y*pi*2) + expr = as_vector([expr, expr]) + expression = as_vector([x, y]) + exact = Function(VectorFunctionSpace(mesh, 'CG', 5)) + exact.interpolate(expr) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V2 = FunctionSpace(mesh, elem.to_ufl()) + res2 = project(function(x), V2) + diff[i - 1] = res2 + else: + V = FunctionSpace(mesh, elem_code, deg) + #res = project(function(x), V) + res = project(expr, V, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'}) + diff[i - 1] = sqrt(assemble(inner((res - exact), (res - exact)) * dx)) + + + #assert np.allclose(res1, res2) + + print("l2 error norms:", diff) + diff = np.array(diff) + conv1 = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv1) + + #print("fuse l2 error norms:", diff) + #diff = np.array(diff) + #conv2 = np.log2(diff[:-1] / diff[1:]) + #print("fuse convergence order:", conv2) + + assert (np.array(conv1) > conv_rate).all() + #assert (np.array(conv2) > conv_rate).all() def test_create_fiat_nd(): cell = polygon(3) @@ -282,6 +335,14 @@ def test_create_fiat_nd(): #nd_fv.to_fiat() +def test_cg3(): + tri = polygon(3) + elem = construct_cg3(tri) + for d in elem.generate(): + print(d.to_quadrature(1)) + breakpoint() + + def test_int_nd(): tri = polygon(3) deg = 2 @@ -294,4 +355,6 @@ def test_int_nd(): dofs = DOFGenerator(xs, S2, S2) int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) dofs = int_ned1.generate() + for d in dofs: + print(d.to_quadrature(1)) breakpoint() From 10e8bbf22d4c71d4e3fa33d6ec362eaa322b2e01 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 17 Nov 2025 12:35:23 +0000 Subject: [PATCH 09/12] fix ordering issue with generation of center dofs --- fuse/dof.py | 10 +++------- test/test_orientations.py | 10 +++++----- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/fuse/dof.py b/fuse/dof.py index 610304d..a6a87fc 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -50,6 +50,8 @@ def permute(self, g): if self.entity: res.entity = self.entity # res.entity = self.entity.orient(g) + if self.orientation is not None: + g = g * self.orientation res.orientation = g return res @@ -96,11 +98,8 @@ def permute(self, g): res = L2Pairing() if self.entity: res.entity = self.entity - #res.entity = self.entity.orient(g) - #if self.entity.oriented: - # breakpoint() if self.orientation is not None: - g = self.orientation * g + g = g * self.orientation res.orientation = g return res @@ -318,7 +317,6 @@ def to_quadrature(self, arg_degree): bvs = np.array(self.cell_defined_on.basis_vectors()) new_bvs = np.array(self.cell_defined_on.orient(self.pairing.orientation).basis_vectors()) basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) - print(basis_change) else: basis_change = np.eye(self.cell_defined_on.get_spatial_dimension()) pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change) @@ -326,8 +324,6 @@ def to_quadrature(self, arg_degree): # need to compute jacobian from attachment. pts = [self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts] immersion = self.target_space.tabulate(wts, self.pairing.entity) - #if basis_change: - # immersion = np.matmul(immersion, basis_change) wts = np.outer(wts, immersion) # pt dict is { pt: (weight, component)} diff --git a/test/test_orientations.py b/test/test_orientations.py index 566252d..574fd78 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -21,7 +21,7 @@ def construct_nd2(tri=None): int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) v_2 = np.array(tri.get_node(tri.ordered_vertices()[2], return_coords = True)) v_1 = np.array(tri.get_node(tri.ordered_vertices()[1], return_coords = True)) - xs = [DOF(L2Pairing(), PolynomialKernel(v_2 - v_1))] + xs = [DOF(L2Pairing(), PolynomialKernel((v_2 - v_1)/2))] #xs = [DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[0]))] # DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[1]))] center_dofs = DOFGenerator(xs, S2, S3) @@ -188,7 +188,6 @@ def interpolate_vs_project(V): exact = Function(FunctionSpace(mesh, 'CG', 5)) expression = as_vector([x, y, z]) f = assemble(interpolate(expression, V)) - breakpoint() expect = project(expression, V) exact.interpolate(expression) return sqrt(assemble(inner((expect - exact), (expect - exact)) * dx)), sqrt(assemble(inner((f - exact), (f - exact)) * dx)), @@ -210,7 +209,7 @@ def test_convergence(elem_gen,elem_code,deg,conv_rate): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, elem_code, deg) - diff_proj[n-1], diff_inte[n-1] = interpolate_vs_project(V) + diff_proj[n-min(scale_range)], diff_inte[n-min(scale_range)] = interpolate_vs_project(V) breakpoint() print("projection l2 error norms:", diff_proj) @@ -305,7 +304,9 @@ def test_create_fiat_nd(): dofs = nd.generate() e = cell.edges()[0] s3 = S3.add_cell(cell) - g = s3.get_member([0,2,1]) + g = s3.get_member([1,2,0]) + print(dofs[-1].to_quadrature(2)) + print(dofs[-1](g).to_quadrature(2)) for d in nd.generate(): print(d.convert_to_fiat(ref_el, deg, (2,)).pt_dict) print(g) @@ -331,7 +332,6 @@ def test_create_fiat_nd(): nd1.to_fiat() nd.to_fiat() - breakpoint() #nd_fv.to_fiat() From 1da3057093a7bb8f9da22f110a12d52861d4be8d Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 1 Dec 2025 15:19:32 +0000 Subject: [PATCH 10/12] test for two forms --- fuse/cells.py | 2 +- fuse/triples.py | 6 +- test/test_orientations.py | 223 ++++++++++++++++++++++---------------- 3 files changed, 135 insertions(+), 96 deletions(-) diff --git a/fuse/cells.py b/fuse/cells.py index 17e6fd4..171b273 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -545,7 +545,7 @@ def get_node(self, node, return_coords=False): if return_coords: top_level_node = self.d_entities_ids(self.graph_dim())[0] if self.dimension == 0: - return [()] + return () return self.attachment(top_level_node, node)() return self.G.nodes.data("point_class")[node] diff --git a/fuse/triples.py b/fuse/triples.py index 617e99c..850c305 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -127,12 +127,14 @@ def to_fiat(self): self.matrices_by_entity = self.make_entity_dense_matrices(ref_el, entity_ids, nodes, poly_set) self.matrices_by_entity = self.dummy_function(ref_el, entity_ids, nodes, poly_set) # TODO remove mat_perms, entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) - self.matrices = mat_perms + if entity_perms is None: + self.matrices = mat_perms + self.reverse_dof_perms() + print(entity_perms) # if not pure_perm: # self.matrices = mat_perms # else: # self.matrices = entity_perms - self.reverse_dof_perms() form_degree = 1 if self.spaces[0].set_shape else 0 # TODO: Change this when Dense case in Firedrake diff --git a/test/test_orientations.py b/test/test_orientations.py index 574fd78..320094e 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -22,8 +22,7 @@ def construct_nd2(tri=None): v_2 = np.array(tri.get_node(tri.ordered_vertices()[2], return_coords = True)) v_1 = np.array(tri.get_node(tri.ordered_vertices()[1], return_coords = True)) xs = [DOF(L2Pairing(), PolynomialKernel((v_2 - v_1)/2))] - #xs = [DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[0]))] -# DOF(L2Pairing(), PolynomialKernel(tri.basis_vectors()[1]))] + center_dofs = DOFGenerator(xs, S2, S3) xs = [immerse(tri, int_ned1, TrHCurl)] tri_dofs = DOFGenerator(xs, C3, S1) @@ -161,7 +160,7 @@ def test_surface_vec(elem_gen,elem_code,deg): assert np.allclose(0, res2) -def interpolate_vs_project(V): +def get_expression(V): mesh = V.mesh() dim = mesh.geometric_dimension if dim == 2: @@ -173,12 +172,10 @@ def interpolate_vs_project(V): if dim == 2: if len(shape) == 0: exact = Function(FunctionSpace(mesh, 'CG', 5)) - #expression = x + y - expression = x + expression = x + y elif len(shape) == 1: exact = Function(VectorFunctionSpace(mesh, 'CG', 5)) expr = x + y - #expr = cos(x*pi*2)*sin(y*pi*2) expression = as_vector([expr, expr]) elif dim == 3: if len(shape) == 0: @@ -187,10 +184,16 @@ def interpolate_vs_project(V): elif len(shape) == 1: exact = Function(FunctionSpace(mesh, 'CG', 5)) expression = as_vector([x, y, z]) + return expression, exact + +def interpolate_vs_project(V, expression, exact): f = assemble(interpolate(expression, V)) - expect = project(expression, V) + #expect = project(expression, V) exact.interpolate(expression) - return sqrt(assemble(inner((expect - exact), (expect - exact)) * dx)), sqrt(assemble(inner((f - exact), (f - exact)) * dx)), + print(exact.dat.data) + #return sqrt(assemble(inner((expect - exact), (expect - exact)) * dx)), + print(f.dat.data) + return 0, sqrt(assemble(inner((f - exact), (f - exact)) * dx)), @@ -199,7 +202,7 @@ def interpolate_vs_project(V): def test_convergence(elem_gen,elem_code,deg,conv_rate): cell = polygon(3) elem = elem_gen(cell) - scale_range = range(2,6) + scale_range = range(0,1) diff_proj = [0 for i in scale_range] diff_inte = [0 for i in scale_range] for n in scale_range: @@ -209,14 +212,18 @@ def test_convergence(elem_gen,elem_code,deg,conv_rate): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, elem_code, deg) - diff_proj[n-min(scale_range)], diff_inte[n-min(scale_range)] = interpolate_vs_project(V) + print(V.cell_node_list) + x,y = SpatialCoordinate(mesh) + expr = cos(x*pi*2)*sin(y*pi*2) + expr = as_vector([expr,expr]) + _, exact = get_expression(V) + diff_proj[n-min(scale_range)], diff_inte[n-min(scale_range)] = interpolate_vs_project(V, expr, exact) - breakpoint() - print("projection l2 error norms:", diff_proj) - diff_proj = np.array(diff_proj) - conv1 = np.log2(diff_proj[:-1] / diff_proj[1:]) - print("convergence order:", conv1) - assert all([c > conv_rate for c in conv1]) + #print("projection l2 error norms:", diff_proj) + #diff_proj = np.array(diff_proj) + #conv1 = np.log2(diff_proj[:-1] / diff_proj[1:]) + #print("convergence order:", conv1) + #assert all([c > conv_rate for c in conv1]) print("interpolation l2 error norms:", diff_inte) diff_inte = np.array(diff_inte) conv1 = np.log2(diff_inte[:-1] / diff_inte[1:]) @@ -232,16 +239,21 @@ def test_convergence(elem_gen,elem_code,deg,conv_rate): def test_interpolation(elem_gen, elem_code, deg): cell = polygon(3) elem = elem_gen(cell) - mesh = UnitSquareMesh(1, 1) + mesh = UnitSquareMesh(2, 2) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, elem_code, deg) - err1, err2 = interpolate_vs_project(V) - assert np.allclose(err1, 0) - assert np.allclose(err2, 0) + expression, _ = get_expression(V) + expect = project(expression, V) + #f = assemble(interpolate(expression, V)) + print(V.cell_node_list) + #print(f.dat.data) + breakpoint() + assert np.allclose(f.dat.data, expect.dat.data) +@pytest.mark.xfail(reason="issues with tets") def test_projection_convergence(elem_gen, elem_code, deg, conv_rate): cell = make_tetrahedron() elem = elem_gen(cell) @@ -283,78 +295,103 @@ def test_projection_convergence(elem_gen, elem_code, deg, conv_rate): assert (np.array(conv1) > conv_rate).all() #assert (np.array(conv2) > conv_rate).all() -def test_create_fiat_nd(): - cell = polygon(3) - nd = construct_nd2(cell) - nd_fv = construct_nd2_for_fiat(cell) - ref_el = cell.to_fiat() - deg = 2 - - from FIAT.nedelec import Nedelec - fiat_elem = Nedelec(ref_el, deg) - print("fiat") - for f in fiat_elem.dual_basis(): - print(f.pt_dict) - - #print("fiat: fuse's version") - #for d in nd_fv.generate(): - # print(d.convert_to_fiat(ref_el, deg).pt_dict) - - print("fuse") - dofs = nd.generate() - e = cell.edges()[0] - s3 = S3.add_cell(cell) - g = s3.get_member([1,2,0]) - print(dofs[-1].to_quadrature(2)) - print(dofs[-1](g).to_quadrature(2)) - for d in nd.generate(): - print(d.convert_to_fiat(ref_el, deg, (2,)).pt_dict) - print(g) - for d in nd.generate(): - print(d(g).convert_to_fiat(ref_el, deg, (2,)).pt_dict) - #nodes = [d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] - #new_nodes = [d(g).convert_to_fiat(ref_el, deg, (2,)) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] - #for i in range(len(new_nodes)): - # print(f"{nodes[i].pt_dict}") - # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") - # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") - # for g in S2.add_cell(cell).members(): - # print(d(g)) - # print(f"{g} {d(g).convert_to_fiat(ref_el, deg).pt_dict}") - - nd1= construct_nd(cell) - print("nd1") - for d in nd1.generate(): - print(d.convert_to_fiat(ref_el, deg).pt_dict) - print(g) - for d in nd1.generate(): - print(d(g).convert_to_fiat(ref_el, deg).pt_dict) - nd1.to_fiat() - - nd.to_fiat() - #nd_fv.to_fiat() - - -def test_cg3(): - tri = polygon(3) - elem = construct_cg3(tri) - for d in elem.generate(): - print(d.to_quadrature(1)) - breakpoint() - -def test_int_nd(): - tri = polygon(3) - deg = 2 - edge = tri.edges()[0] - x = sp.Symbol("x") - y = sp.Symbol("y") - - xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] +@pytest.mark.parametrize("elem_gen,elem_code,deg", [ + # (create_cg2_tri, "CG", 2), + # (create_cg1, "CG", 1), + # (create_dg1, "DG", 1), + # (construct_cg3, "CG", 3), + (construct_nd, "N1curl", 1), + ]) +def test_two_form(elem_gen, elem_code, deg): + cell = polygon(3) + elem = elem_gen(cell) + mesh = UnitSquareMesh(2, 2) - dofs = DOFGenerator(xs, S2, S2) - int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - dofs = int_ned1.generate() - for d in dofs: - print(d.to_quadrature(1)) - breakpoint() + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V2 = FunctionSpace(mesh, construct_nd2(cell).to_ufl()) + else: + V2 = FunctionSpace(mesh, "N1Curl", 2) + v = TestFunction(V) + u = TrialFunction(V2) + assemble(inner(u, v)*dx) + +#def test_create_fiat_nd(): +# cell = polygon(3) +# nd = construct_nd2(cell) +# nd_fv = construct_nd2_for_fiat(cell) +# ref_el = cell.to_fiat() +# deg = 2 +# +# from FIAT.nedelec import Nedelec +# fiat_elem = Nedelec(ref_el, deg) +# print("fiat") +# for f in fiat_elem.dual_basis(): +# print(f.pt_dict) +# +# #print("fiat: fuse's version") +# #for d in nd_fv.generate(): +# # print(d.convert_to_fiat(ref_el, deg).pt_dict) +# +# print("fuse") +# dofs = nd.generate() +# e = cell.edges()[0] +# s3 = S3.add_cell(cell) +# g = s3.get_member([1,2,0]) +# print(dofs[-1].to_quadrature(2)) +# print(dofs[-1](g).to_quadrature(2)) +# for d in nd.generate(): +# print(d.convert_to_fiat(ref_el, deg, (2,)).pt_dict) +# print(g) +# for d in nd.generate(): +# print(d(g).convert_to_fiat(ref_el, deg, (2,)).pt_dict) +# #nodes = [d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] +# #new_nodes = [d(g).convert_to_fiat(ref_el, deg, (2,)) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] +# #for i in range(len(new_nodes)): +# # print(f"{nodes[i].pt_dict}") +# # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") +# # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") +# # for g in S2.add_cell(cell).members(): +# # print(d(g)) +# # print(f"{g} {d(g).convert_to_fiat(ref_el, deg).pt_dict}") +# +# nd1= construct_nd(cell) +# print("nd1") +# for d in nd1.generate(): +# print(d.convert_to_fiat(ref_el, deg).pt_dict) +# print(g) +# for d in nd1.generate(): +# print(d(g).convert_to_fiat(ref_el, deg).pt_dict) +# nd1.to_fiat() +# +# nd.to_fiat() +# #nd_fv.to_fiat() +# +# +#def test_cg3(): +# tri = polygon(3) +# elem = construct_cg3(tri) +# for d in elem.generate(): +# print(d.to_quadrature(1)) +# breakpoint() +# +# +#def test_int_nd(): +# tri = polygon(3) +# deg = 2 +# edge = tri.edges()[0] +# x = sp.Symbol("x") +# y = sp.Symbol("y") +# +# xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] +# +# dofs = DOFGenerator(xs, S2, S2) +# int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) +# dofs = int_ned1.generate() +# for d in dofs: +# print(d.to_quadrature(1)) +# breakpoint() From b9b8c06a7ca9f2eb15e600d5647170af98c9fdc3 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 6 Jan 2026 11:45:05 +0000 Subject: [PATCH 11/12] add 3d polynomial example --- fuse/dof.py | 2 +- test/test_polynomial_space.py | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/fuse/dof.py b/fuse/dof.py index a6a87fc..48b857d 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -318,7 +318,7 @@ def to_quadrature(self, arg_degree): new_bvs = np.array(self.cell_defined_on.orient(self.pairing.orientation).basis_vectors()) basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) else: - basis_change = np.eye(self.cell_defined_on.get_spatial_dimension()) + basis_change = np.eye(self.cell_defined_on.get_spatial_dimension() + 1) pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change) if self.immersed: # need to compute jacobian from attachment. diff --git a/test/test_polynomial_space.py b/test/test_polynomial_space.py index cb1c769..0da6486 100644 --- a/test/test_polynomial_space.py +++ b/test/test_polynomial_space.py @@ -152,6 +152,41 @@ def flatten(coeffs): nc = np.reshape(coeffs, newshape) return nc +@pytest.mark.parametrize("deg", [1, 2, 3, 4]) +def test_3d_nd_construction(deg): + cell = make_tetrahedron() + ref_el = cell.to_fiat() + sd = ref_el.get_spatial_dimension() + x = sp.Symbol("x") + y = sp.Symbol("y") + z = sp.Symbol("z") + M1 = sp.Matrix([[0, z, -y]]) + M2 = sp.Matrix([[z, 0, -x]]) + M3 = sp.Matrix([[y, -x, 0]]) + + vec_Pd = PolynomialSpace(deg - 1, set_shape=True) + Pd = PolynomialSpace(deg - 1) + composite = vec_Pd + (Pd.restrict(deg - 2, deg - 1))*M1 + (Pd.restrict(deg - 2, deg - 1))*M2 + (Pd.restrict(deg - 2, deg - 1))*M3 + + assert composite.set_shape + assert isinstance(composite, ConstructedPolynomialSpace) + on_set = composite.to_ON_polynomial_set(ref_el) + + from FIAT.nedelec import NedelecSpace3D + fiat_nd_space = NedelecSpace3D(ref_el, deg) + + Q = create_quadrature(ref_el, 2*(deg + 1)) + Qpts, _ = Q.get_points(), Q.get_weights() + fiat_vals = fiat_nd_space.tabulate(Qpts)[(0,) * sd] + my_vals = on_set.tabulate(Qpts)[(0,) * sd] + fiat_vals = flatten(fiat_vals) + my_vals = flatten(my_vals) + + (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) @pytest.mark.parametrize("deg", [1, 2, 3, 4]) def test_3d_rt_construction(deg): From ac11f64d780d5196087bb75e8b679e9140dda3e7 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 6 Jan 2026 20:12:14 +0000 Subject: [PATCH 12/12] some tidying up, fiddling with tests --- fuse/.triples.py.swp | Bin 0 -> 40960 bytes fuse/dof.py | 10 +++-- fuse/traces.py | 1 + fuse/triples.py | 17 ++++---- test/test_2d_examples_docs.py | 27 ++++++------ test/test_3d_examples_docs.py | 75 +++++++++++++++++++++++----------- test/test_convert_to_fiat.py | 9 ++-- test/test_orientations.py | 33 ++++++++------- 8 files changed, 102 insertions(+), 70 deletions(-) create mode 100644 fuse/.triples.py.swp diff --git a/fuse/.triples.py.swp b/fuse/.triples.py.swp new file mode 100644 index 0000000000000000000000000000000000000000..c6ec9c8735dee8f41a5a02a5ecf92d19c46cc38c GIT binary patch literal 40960 zcmeI53y@@2dEZ+$7$Jd;kzK^tflCi~d(doaXJiZvS}bW-D{Vm1qFqbEv!kY``%X{q z_DkJ8yEER5aA6~`olv%kjW85(j1!lg3dPHIxqwQ51P2i#gdB)tY|6n-3PPNUow%GJ zLw^7BxcA(?w`XQoS`sIHYyQ)>?>*1&ob#RYJESj3d|_*Hitg`*HXAT`mKk&(qxVSGdpBiT78!=gaQ%p^5iTchB#4pP!w0 z|E=!%oclcDD(v%jk7pk<3d|@lqri*;GYZToFr&bX0y7HCC@`bIi~=(X{69&7db3pe z4r=ys3BuU_z5M?n9&Uj{;Mw3)-(D*HBk+1~0=y7h4IY1Hsq_W#yWpRKcYrs83!nz> z1FPT|cqRB&@Gam^o>40O6u1$*47?Ov34Z(PQt79`S@0TA22TY~0sr}`Qt7kcH^5o& z8t@GS|3|<_z+1r?a11PfYr!?(_fa0a7TgPN2QL9H1kVItLhJL8 zFZtKQc)b>Pe1JyW=~T!4+IqV_c7dAldgPxdk$H42Y2@+duooS}P2-kM+>LwViV0!@ zxQLtU9s4kdr})Fz#0`oPcSiBDiK_Ryd>B)BQ?^z-u8juual78BHQVGkYWI7U`GvGX z?R8sMLN(&i%AwWjNO3mmgCOy_=+vD@?~G3BX*AmCZ+4o|8h3-Z(Oz#i;$}4JcjNOL z@i2~5tY|RouhrK&+a(oCHB&*()jON9O*&UaSM&xniH8eOqu=Ro4o55F_Sqj@U7<~? zs!olR+P%10>-VZjZB|qIDl?k`4b*9mAn#_k(i>Fk!(n~fY-^2v zf7l!yJ~iBo=ebHsNLQZ3N%AR;P08xxalNr&nzX32?dE*t&_Z-5Od84;1x+dz9X=c_ zWgGa3)xN8mLxX1nV zCP2M$w!E`ba&dicqElQrR@bO;3Helf8a{q`XTikv9hYXRwR`J*(}Nb&Y73EirtL7n zKy|&?^o?{tHDi@77`bjZjW^n(s6FZ(7&CI~oz6A|9S-BhILX|_G)$r(8#=efMh!ky z?QM2z+>a{r(RzQVGU{>E?%Cme^v+|qT5i@KuG$7s?UGoNS!}ms#<;0S7&Rej@ZMIm z=#553QkPn-&CssZrVagJJWA4txY|vI<)+v{bSyhPJ?Aqs%V^FGTaqgQOo8JzZQadb*^i<$KF%tzE}1Yw9{>dn;AXU`TyEH=1ye$v$uXSVtsPZ*r*;T3 zu4k2kyQ4D=tRHW8_!2jC`in{It7LM!xO`PbF(qF;Y!NbcKUX&J$9X*)r?8Y zlHE)nV$FVGST2`MyyZcCSnn?TN-amXdVbC1CvRSsl`Nwyub{bOa?-pqiu&u(IycRJ zx0*@GX1MI56U;xUoalG9d;KoU!jo$7g~;^3zU_Yt(Q&ceL7%}%$1x%7F!5zng_ffu z{6ErH-`ujTe?;9~#nh0G>4IZZ&hZAN)7-|ous$4{0;~M7?bN2sqJkM0en2=bBPx~F zhy8A}zS(Px`~416W_Qr1|5z63zKj}ee4X7FX0^dK-7M-2{3HCMd+lCpfs_YJ$$&@N^CX94e1Jj(Zz@z~Xx zVpccD?UaOuP&LGr(Rf%hOf(4D8um8_3AG3Duv>!&!}gXUn6GA8y%3puEw}Z;1__cG z55;YKflldL@Ro7wUmcA*AI?Adv?_k zv^KKrRI*QyxuaLeYlEVwQzzv#19~CK4|2&SUAl2hWf>v;RvT6-`M-o5{VPE7{~RB> z-i^%P0TK8^mFwLq+w7WAU`Bx%1!feOQD8=a83kq(`2U6ih^EM8 z#)xpd+l_~#xQWa?s*LMHDd(+>%JLp@4(+CWXqL|lQLo>O=Svc7+cxHYBON|ua+8Hbpjqyr3C`cLQTkygZYVsKi`*NVktVu{TDOG$-FJT8+mxHQXv0JrL{+n)E*ZyMP9Zlk)Rgn&@4U6{+**Pk5&5QDPOFz+UO1tqEw)1r7`L> zxf7Yvt(Lh%i=HXmXq9NqORV+VbXY#m92W6!UQ2l1cd#sXlK>btp+0 zC_+lxsroOtZaq0AjZKHat zr8XXG9`KpE($OQldZwu=RsCSGIj6+Nt~jn~O!;AC##frA=$uI@(c;_g6xbyLqku}| zXqpzZ^5TAFE6`A?gB8@L^5U&So_fn4oCf} ziCJy4>3-ZEZ>PP~jZcvr9#d2ksWUdEb5TsB`-<;8E~ohwlP{{vp>4)hBmdvVO!8r8 z|M&j=pGWThb?|dQegHAJ30wny54ry_@Btt_z`fue@O5PVuYiZad%&B(4d4%v`5y$Y z1+w)&ifrEkSAZ`f+kXVS6Z{y6!4H8$;M>4oW553x_$YWYXn^N~?*xB<{r(~FHgF0Y z1PkB_@G!RfPlC6AH-Hwn8${p>*zErRd>Fh1Y=WD?Gr*(R?%xis1K$IF2HX7#xB+|) zTm8QRKMP`T8@Lrb4?G(@hK>G%;6d<4@H+6l;JM&0q3`E`==`gq_1;(5^}_r57HyCI zu3Gr8)u;#dIt5v8w_F77Hiw}(&-X+Pfz&P?-GDJNs}&#Sx_)%j8|CS}Rsq$h_j(Kn z%`TPX{KQ|;#P;IzJgy@dvQRP`EdmB-(U@u)s#Mp}TH81(Ah_;7g~cL&u$uT>5ikmF;r zv0?=5b}^l#_6jBI0+dhp%GERI=qs6^Y2_3)^GXxs%7N3TvHx7qivcq&t+{g`V>!UZ z;e5(6l!&f;LgMm{p_rwTd_fedVjK#Ls(;dPfr*HNJ>WfZaTw00+1flXRLK{Crn9#G zl7FnvxRbMLR0lbQDTE9o`TKn7*UP6*(-58Juzr5!f+^3=YIGs56!z6IC7<7k8tq{N zdqibvfNeyabOkHU2D29PLVt6-K_VA+c3x4=SFYH;rXuExcR7fgr(f87=zFy=D4hVVcTPpG?Z4pp&5>AajV_i_kheMkr>mJ zy5~2#zbn3fvn_MGxYX&>m;XAgip|sNzGyr!W)o-`Epj&V5Wh_DkB!ObrP}XL3fRipwxx zmNZs=yeq5yq>bc!vd!9i|hAOmq5|X;zX73tdtGU2|Fb zDB@&1dWJE!rSsn$Wm1&~%v)&t7%}@v+WS&4u)6(o6GWSyQI88|&H76uq#ARXJmw{mcCwK`y8U=@=NkuWtgM!gu# zv@u#~uU7F1kq4yAwZ?PRm(??Xz)>@h9m8Z@FM#i?AblARO`Ec$Kp3I-(rh;o;=1jg zIFt#w>dWG)nJ;xRz*iPmQ8tQW}-Y z2DH4u1R$FWJ4Hb}treiXjH_8n46-$G(|D?=h>$H>AWisLg4E2wO*tXm=gvE7dT(7x zH0!v2pEj#&d>`Aa7mQc3#Ao~gDQv)lC?6qhe>7?mXJp?F!k8{iDCa^qNgLh-WN}m@ zm+bSPi}S)%Y|=$Jt%DsGK3qxw`<$Z5{*RblIwzSF`@i46{|V&$Ujsw%Ebu7u{X<|A zyapTq{}tK(A@Bxp2l#$)75D-&{(HdBgFcY$|7P&L;BjRBhk8?}zVU#pLpqFnpRZ$o;WfKhKcNM-t9=Msu+XjClo^zG1Ud&@ukb^Y1qWw=PycF^|~NCo6+tPHc9xTM$qa%|A5(n|iRVS=WFn`@fp0{*8UxH6ijlBvxI zvT6xCAxd$Uchq&H6TO0h8Xly4EU(5?$O3B4rjd(DQD#kl6o7O$nAy1MK zo_pDMCZC|aeaR+DP4j0WkNTZ+afR`JGX^ZW(rD1f$Yq`MOFj*Kj)#FlSP#5+U=?); zlal{hDaX=K2bDVC9&eNSid7rceMqc zfo9`Pr83$Wz*QEagNs4nn+PoN2~zJ7lITHLlV-Y@a7nlTQdwc(82xA>4yN;7S~3Fr zb%2?cZF!mbN+O@G{mOMVb@lE-etb#&TKVU7t0$yBl}SOu(r9xar%`q)!Dfs+2+Na? zbiWrgZ&n;rvbv-tmM9;rLG-vdFxuhjsD6$maI50AoGC>W5BBa~r0Yq&sBanCZr8Wk ztclD%taG>>uD6b^({F%%W-&$-qR=;51aoB+RxD%ZyY{UT`JDQE8GgSpen`_(u67j>Wgj)_4+WW za!|hFG8H$tQq#2yRix3~n<&MZ(pv)8^T|x6Yb;eHPejGK2(6`}=Cj?8=29_r*;a#T ztZ364f9X_d1DFh4sGJhP2l2hhGz^I)^E34+(1R+PSxLyH*z6dXT(2-^nsrrG9rk+c z{eix+zIC`BSVfeWF%fAzYcB(<3)9Wn*?5M{9q^6A2d-~Uj3fe!Q!>`L%NmJHk2N*0 zt0d2Tw6MEA%LkPAjOU07>2%6sCs^kS{sX*^n+cZ zo&e+^4TcFkuP!sOJJW~rsa+1X>C?W9+{EpXO7{i9xXv0^|BT;%!P4;Z@`?S{Pt|9q<4w*;|m%tc| zKm`6Anf)W+XTV8t7kD9%{r(Z~>)0l1*w*fiVXn zXW|6IVmh$8WCNLbrYW5+&$IGhVf-q$p(Zf8BRj@N;5w*o0-Ze1nkZ}B6ZEc zL~G$32lMllaw^C{*7kxdU}vH&IuFg%Nj%2-wl**$RvMioe9pbg7oZf@&EJa(#+lHi zzGoG0NHahWDnVhdO6IfJR~>z{u!_JO2s)=0&IsJvY8ZJ%3xr8;he1#@m#wq&7Mydw~n^`;mC{>_t==ZxBbgFmCRMN zu$5(yE_t#~eUS{xDdbA199m8ML*h)Z<|IR(@x*7)n$ubPPLfJN1Cr+jD9zap1#am{tkL28QWIX7GJwQr zaTyI5`_eqCTApUXpu&)#Q4}1jvJcxr9CCt^Qr2arlBpcx;)b5lkz=wSa4gypGO zbG^yzo33g6YDQgITsX9v_5+`a|MUfULJVPr)hipfWdfHC$|$3BQ{aE-^9^~#H;-lb z+@Czwp~je;(M5K5Qq|=&3;D|EvzRye5W6?B=TWZpCb*PO!>rr);#WSfK6lA9u&swN zpru*f*Fke>KSOrj8gxrOC4FjTvKJ4R{C^!n?Ddjmk^lYq`=3Sb|7YOGfPDX#!MB6Y zA@~0(cq3Q=Is@=A^Z@S!I{WWl@N#f9_yBr=p93}U-_ZeRum5j?*MV!mAEN_!KX@ZJ z1LXIA4ft>90e%O(6Wj%)7kDqafH#1@1%4kLz!{)3`aX`_FFXF{koE6!^8PO)+jqcU zAkV)8{2lOJx<`)R0Y3=-BlUd`P`#H__x)cXi6`&p(mg6!@qCg8dbS2@H=BqoZwtT8 zYvH67FbT2x7=(@_liBL~oDWR2R+UIP??RlqfX|2Vy%Dggn@aiE51+hcIkT6=_#A3T zG&EUdy!yj*{JlFO?ED702gXyqmq`qjxf@-#7W;@PncX>#LQGMDTz)FM2DwNYBO=Lh zQv*`6lT76DGmALeNHUO?S7y33q;Z;xlImG>8E)yB84^{^oIc|?la)|SVs0W&SqBd) z6O!d*`hVl-dcrJF?4~OOO!kF(qS{51=*xjR9N?4Eyfj&B|J}|*&IBVJk0q_?Rn^vF zss~8udqcpPzv<)+^m-Y_kr!vO70d|+45z4)%F|L}sefe>)f^XkviKy!&nLO}$mnJG zSjN44a*UV45=`F1k9(X>h=ZIpox_j{I)ojBo19?Cxlc5swXW}xFiB|_vbvk=YU*Vg z!J493sM<`uf63I;NJo)!1KRbvoc1Z{do&0FhgldWUDb1e{7S@AMk}!-IBd)}GRlxZ zYc;jB=H5W;fmf4mRCvIWZO~EIs3BbN{Og|SrtBdl5)aP8}t7g1W%`CY~&kp8M zYQ@m7C#SAVD;h2;&GGgiKCENXrdBtpNU`RPzOw2(=^R~)>*E-|R~wTmJ3DKkN<8S* zt__QuFh@WtJ7vkcyMXUWoh%ZaV4u&{y-;77)B)t^whzHt87otx&w>(xCQ+~Sgsg&D zD_!bB?zE{4O|ARybiG?zn>o$AiGH%)U|Vieno6ggWJu`1WVF5|l|t8?EhfiVx&e;w z2fuF}Z4SJulF)XxVX0%=)C^9yL>E^nnUwf9g*CEeWzy7RYChG}qUFPUsv$KZy_5S* z#5})TFdH2LlS@*jg2-$Nn4*!`k1k@5tjtNwcrf+OM)fsz+i7D9EYtWjv9qufwfe!~ zjvOAz#!GJ}(ry`JC`enhU=#(loG>i;E@eAYcI;XA`h+a|*W##g)gV4`*NOZNTiMKv zl1IG3ObSm(ud~lCl+;6|^H;W)1KtsDSq%=idvS10F@rmoNS)a25Cv zbsT_$;A)^cKOMNscNTr36c?v3s~fiz4tDZ|Ar04co6}79N0-iN^y$g~I-9YNr!41u zhEc=&G>9LOEMXW-vYGg;Z{2>AGwbE0$cuC$43p8M$QUyr67BtNrK@sAXBW$6%2vbV z(o#{2J!`s~mLfAbJ*m{WRI2!yGAPW&5~C=IT{Jb5?dWM*kj%<#(&c-dW{qC81%L65=ak z%Ba&)3Us*tXO^Iu{+(OFd*VW=^f>i`I`Te$Mr*|%u4`5g)(sTbcm3Q3W=T{>I*tP0 z1ml98?m7#%qNod`Z&KEAaZ4TRsD4L(d^|nX+ik2krG3gcrRwK3B*rID%Q4QlzU$8sav?VW7X1d0127-8Npsh@@1@37fP_gFN4Z{yk z$c+!Zm^w7(U8DbE0ng7yaEvs*JTc3{qDq9%e9m+fmMUk5obb?+uek<@6~4F$MAa~E zc%N$2r;*6}?#Z2R)p$GErd7-di`uGhYP}QbsF5v<`?8DKc(PC~hE;I+ZI4E4aV|e+ zg^_iyAzL}S5S?RFq*oAT6+TmHl^@2FN&D|X6={PhN+2ItpC$qnA+)*8dxyoY6PChg zMfWDY(n&TZXo`_13O`%~<>xjy-`(@bK1i0usijF!*Ny&xqRO=Uy@b7 zsagn}6JDfT7(-ZR9m0qdKcA2HA8Nm2+ zx_MtJ8AX;tg?u_5sF$RwsmY3T12n+924H422d@>axz@{fn&TafnNY!a2$LexB|!r;O(FbUIo62 z4&XO|&ii{1^uTk#6`%yZj9%bT@EPzy5Q8RYfZM^V!PVeP=nOstbT+^PKz;%@f_YE@ z&jMdVkD%WLcm#YJycfI)+y_>{Eg%A4M#u0;@UOw!fc6)xfg|7`xB}=rfnNqga2wD$ z0*`?Yf(OBm0PQ(=2Dlo06+Ohi1s?%_5BxN!gV%u@!Ij`?;7`y?{1$i!+zVa_4ucl~ z=`6m2p5nK`kAt(|47dw?H~0p6i?4(K1RepO2JZzrZ?Fj(;CAq8a25Cx?e+=qE^rRK z2AlwD&mRDn0zddo_;fF*h*!70VfC|?^C5{yioPi zGiiMeBRMlZ7LJzT%^ep0(M_?lu=26bgAZ34i+kz0Pf#Jc}e+6IE=@hk7J66 zrP(Y_vRa9Bn50^XTyQhua%bHkmYzQLBWkE3~J%Xr4SK zDUwwild}1JPk$hazzwst0~&R#zd=tJwXEm38s$p6&D&wEC7)?*5Q&;r`E^36d1*UagbyU~mO3L0}`4cCMCyPkxGi*lN z@i;hN!WS|~QtQoZX{F{QfOTIF3uWrpTB6?W&UP^7_=z9cv#%GEfQv-|MVRDA`V_W9D ztv2*8dZ#0wUDNc3=19$SikU0;S;AXuOAFDh`WdvehZl!q0uREus`rAQ0G!BHeod~3 zZa0oCuc&X=Kbhr?9w*gwjMe>^*(tM_s@PS+N?B8aS@GC?kY#&VS@=XyA!o0n@xy0n)2r|3 z%Di=p45la<_N8kWO}9O!uUEQRlaFyf4qfnqbXO`5bSErCDL=DCvb@lo&iABv7NrHH zQie=M5udfSPSb*Z6Kz7p%|V%ib%~i6t?N>zCONxNrd_jXq>ChzM6#gDuaMN?Ttec` IS conv_rate).all() -@pytest.mark.parametrize("elem_gen,elem_code,deg", [ +@pytest.mark.parametrize("elem_gen,elem_gen2,elem_code,deg,deg2", [ # (create_cg2_tri, "CG", 2), - # (create_cg1, "CG", 1), + (create_cg1, create_cg2_tri, "CG", 1, 2), # (create_dg1, "DG", 1), # (construct_cg3, "CG", 3), - (construct_nd, "N1curl", 1), + (construct_nd, construct_nd2, "N1curl", 1, 2), ]) -def test_two_form(elem_gen, elem_code, deg): +def test_two_form(elem_gen, elem_gen2, elem_code, deg, deg2): cell = polygon(3) elem = elem_gen(cell) - mesh = UnitSquareMesh(2, 2) + mesh = UnitSquareMesh(1, 1) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, elem_code, deg) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V2 = FunctionSpace(mesh, construct_nd2(cell).to_ufl()) + V2 = FunctionSpace(mesh, elem_gen2(cell).to_ufl()) else: - V2 = FunctionSpace(mesh, "N1Curl", 2) + V2 = FunctionSpace(mesh, elem_code, deg2) v = TestFunction(V) u = TrialFunction(V2) - assemble(inner(u, v)*dx) + res =assemble(inner(u, v)*dx).M.values + for row in res: + print(row) #def test_create_fiat_nd(): # cell = polygon(3)