diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index b33e9ac..98010e5 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -62,7 +62,7 @@ jobs: with: repository: firedrakeproject/firedrake path: firedrake-repo - ref: refs/heads/connorjward/pyop3-develop + ref: refs/heads/indiamai/fuse_mesh_cell - name: Validate single source of truth run: ./firedrake-repo/scripts/check-config @@ -116,6 +116,7 @@ jobs: : # Fix for petsc4py+slepc4py build echo 'setuptools<81' > constraints.txt + echo 'Cython<3.2.6' >> constraints.txt export PIP_CONSTRAINT=constraints.txt : # Install build dependencies diff --git a/constraints.txt b/constraints.txt new file mode 100644 index 0000000..3f08c3a --- /dev/null +++ b/constraints.txt @@ -0,0 +1 @@ +Cython=3.2.5 diff --git a/fuse/cells.py b/fuse/cells.py index 57d5dfd..4d9b42c 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -289,7 +289,7 @@ class Point(): id_iter = itertools.count() def __init__(self, d, edges=[], vertex_num=None, oriented=False, group=None, edge_orientations=None, cell_id=None): - if not cell_id: + if cell_id is None: cell_id = next(self.id_iter) self.id = cell_id self.dimension = d @@ -299,7 +299,7 @@ def __init__(self, d, edges=[], vertex_num=None, oriented=False, group=None, edg if d == 0: assert (edges == []) - if vertex_num: + if vertex_num and vertex_num > 1: edges = self.compute_attachments(vertex_num, edges, edge_orientations) self.oriented = oriented @@ -891,8 +891,10 @@ def to_ufl(self, name=None): def _to_dict(self): # think this is probably missing stuff + # d, edges=[], vertex_num=None, oriented=False, group=None, edge_orientations=None, cell_id=None): o_dict = {"dim": self.dimension, "edges": [c for c in self.connections], + "vertex_num": len(self.vertices()), "oriented": self.oriented, "id": self.id} return o_dict diff --git a/fuse/dof.py b/fuse/dof.py index 6d8ff2c..27d1847 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -260,14 +260,14 @@ def evaluate(self, Qpts, bary_pts, Qwts, basis_change, immersed, dim, value_shap return Qpts, np.array(wts).astype(np.float64), comps def _to_dict(self): - o_dict = {"fn": self.fn} + o_dict = {"fn": self.fn, "syms": self.syms} return o_dict def dict_id(self): return "BarycentricPolynomialKernel" def _from_dict(obj_dict): - return BarycentricPolynomialKernel(obj_dict["fn"]) + return BarycentricPolynomialKernel(obj_dict["fn"], symbols=obj_dict["syms"]) class PolynomialKernel(BaseKernel): diff --git a/fuse/groups.py b/fuse/groups.py index 0ada9d6..5bd9ad7 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -297,6 +297,16 @@ def __repr__(self): return self.name + str(self.size()) return "GS" + str(self.size()) + def _to_dict(self): + return {"members": [m.perm.array_form for m in self._members]} + + def dict_id(self): + return "PermutationSet" + + def _from_dict(o_dict): + perm_set = [Permutation(m) for m in o_dict["members"]] + return PermutationSetRepresentation(perm_set) + class GroupRepresentation(PermutationSetRepresentation): """ diff --git a/fuse/serialisation.py b/fuse/serialisation.py index c41e642..fcd8e66 100644 --- a/fuse/serialisation.py +++ b/fuse/serialisation.py @@ -31,6 +31,7 @@ def __init__(self): "Edge": Edge, "Triple": ElementTriple, "Group": GroupRepresentation, + "PermutationSet": PermutationSetRepresentation, "SobolevSpace": ElementSobolevSpace, "InterpolationSpace": InterpolationSpace, "PolynomialSpace": PolynomialSpace, @@ -40,8 +41,11 @@ def __init__(self): "DOFGen": DOFGenerator, "Delta": DeltaPairing, "L2Inner": L2Pairing, + "BarycentricPolynomialKernel": BarycentricPolynomialKernel, + "VectorKernel": VectorKernel, "PolynomialKernel": PolynomialKernel, "PointKernel": PointKernel, + "ComponentKernel": ComponentKernel, "Trace": Trace } @@ -65,6 +69,10 @@ def encode_traverse(self, obj, path=[]): res_array[i] = dfs_res return res_array + # Some sympy objects are not hashable, so we must accept we may serialise them twice. + if isinstance(obj, sp.core.containers.Tuple) or isinstance(obj, sp.Expr) or isinstance(obj, sp.Matrix) or isinstance(obj, sp.Poly): + return "Sympy/" + sp.srepr(obj) + if obj in self.seen_objs.keys(): return self.seen_objs[obj]["id"] @@ -75,9 +83,6 @@ def encode_traverse(self, obj, path=[]): self.store_obj(obj, obj.dict_id(), obj_id, obj_dict, path) return obj.dict_id() + "/" + str(obj_id) - if isinstance(obj, sp.core.containers.Tuple) or isinstance(obj, sp.Expr): - return "Sympy/" + sp.srepr(obj) - return obj def get_id(self, obj): diff --git a/test/test_serialisation.py b/test/test_serialisation.py index 56db668..0974522 100644 --- a/test/test_serialisation.py +++ b/test/test_serialisation.py @@ -1,7 +1,9 @@ from fuse import * +from firedrake import * from fuse.serialisation import ElementSerialiser from test_convert_to_fiat import create_cg1 -# from test_2d_examples_docs import construct_nd +from test_orientations import interpolate_vs_project, get_expression +import pytest import numpy as np vert = Point(0) @@ -60,13 +62,36 @@ def test_cg_examples(): assert any([np.allclose(dof_val, dof_val2) for dof_val2 in dofs]) -# def test_ned(): -# cell = polygon(3) -# triple = construct_nd(cell) -# converter = ElementSerialiser() -# encoded = converter.encode(triple) +cg_params = [(0, 0, deg, deg + 0.75) for deg in list(range(1, 3))] + [(1, 0, deg, deg + 0.75) for deg in list(range(1, 3))] +nd_params = [(0, 1, deg, deg - 0.2) for deg in list(range(1, 3))] +rt_params = [(0, 2, deg, deg - 0.2) for deg in list(range(1, 3))] +dg_params = [(0, 3, deg, deg + 0.75) for deg in list(range(0, 3))] + [(1, 3, deg, deg + 0.75) for deg in list(range(0, 3))] +nd2_params = [(1, 1, deg, deg + 0.75) for deg in list(range(1, 3))] +bdm_params = [(1, 2, deg, deg + 0.75) for deg in list(range(1, 3))] -# decoded = converter.decode(encoded) -# for d in decoded.generate(): -# dof_val = d.eval(phi_2) -# assert any([np.allclose(dof_val, dof_val2) for dof_val2 in dofs]) + +@pytest.mark.parametrize("col,k,deg,conv_rate", cg_params + nd_params + rt_params + dg_params + nd2_params + bdm_params) +def test_post_serialisation_convergence(col, k, deg, conv_rate): + "Tests that appropriate convergence is still achieved after serialisation." + elem = periodic_table(col, 2, k, deg) + converter = ElementSerialiser() + encoded = converter.encode(elem) + elem_decoded = converter.decode(encoded) + scale_range = range(3, 6) + diff_inte = [0 for i in scale_range] + for n in scale_range: + mesh = UnitSquareMesh(2**n, 2**n, use_fuse=True) + + V = FunctionSpace(mesh, elem_decoded.to_ufl()) + x, y = SpatialCoordinate(mesh) + expr = cos(x*pi*2)*sin(y*pi*2) + if len(elem.get_value_shape()) > 0: + expr = as_vector([expr, expr]) + _, exact = get_expression(V) + _, diff_inte[n-min(scale_range)] = interpolate_vs_project(V, expr, exact) + + print("interpolation l2 error norms:", diff_inte) + diff_inte = np.array(diff_inte) + conv = np.log2(diff_inte[:-1] / diff_inte[1:]) + print("convergence order:", conv) + assert all([c > conv_rate for c in conv])