Test loaded elements function correctly

.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

constraints.txt

@@ -0,0 +1 @@
+Cython=3.2.5

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

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):

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):
     """

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):

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])