tidy up

Author: indiamai

fuse/__init__.py

@@ -1,5 +1,5 @@
 
-from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex, TensorProductPoint
+from fuse.cells import Point, Edge, polygon, line, make_tetrahedron, constructCellComplex, TensorProductPoint
 from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, diff_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group
 from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, VectorKernel, PolynomialKernel, ComponentKernel
 from fuse.triples import ElementTriple, DOFGenerator, immerse

fuse/cells.py

@@ -926,7 +926,9 @@ def __call__(self, *x):
                     if len(attach_comp.atoms(sp.Symbol)) <= len(x):
                         res.append(sympy_to_numpy(attach_comp, syms, x))
                     else:
-                        res.append(attach_comp.subs({syms[i]: x[i] for i in range(len(x))}))
+                        res_val = attach_comp.subs({syms[i]: x[i] for i in range(len(x))})
+                        res.append(res_val)
+
                 return tuple(res)
             return sympy_to_numpy(self.attachment, syms, x)
         return x
@@ -1046,7 +1048,7 @@ def construct_fuse_rep(self):
                 # of all combinations of point, take those where at least one changes and at least one is the same
                 if any(a[i] == b[i] for i in range(len(a))) and any(a[i] != b[i] for i in range(len(sub_cells))):
                     # ensure if they change, that edge exists in the existing topology
-                    if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i].topology[1].values()) for i in range(len(sub_cells))]):
+                    if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i]._topology[1].values()) for i in range(len(sub_cells))]):
                         edges.append((a, b))
         # hasse level 1
         edge_cls1 = {e: None for e in edges}

fuse/spaces/polynomial_spaces.py

@@ -157,7 +157,7 @@ def __init__(self, weights, spaces):
         self.weights = weights
         self.spaces = spaces
 
-        weight_degrees = [0 if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)) else max_deg_sp_mat(w) for w in self.weights]
+        weight_degrees = [0 if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)) else max_deg_sp_expr(w) for w in self.weights]
 
         maxdegree = max([space.maxdegree + w_deg for space, w_deg in zip(spaces, weight_degrees)])
         mindegree = min([space.mindegree + w_deg for space, w_deg in zip(spaces, weight_degrees)])

fuse/utils.py

@@ -29,7 +29,7 @@ def sympy_to_numpy(array, symbols, values):
     """
     substituted = array.subs({symbols[i]: values[i] for i in range(len(values))})
 
-    if len(array.atoms(sp.Symbol)) == len(values) and all(not isinstance(v, sp.Expr) for v in values):
+    if len(array.atoms(sp.Symbol)) <= len(values) and all(not isinstance(v, sp.Expr) for v in values):
         nparray = np.array(substituted).astype(np.float64)
 
         if len(nparray.shape) > 1:

test/test_convert_to_fiat.py

@@ -563,13 +563,27 @@ def test_quad(elem_gen):
     assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9)
 
 
-# # @pytest.mark.xfail(reason="Issue with quad cell")
-# def test_non_tensor_quad():
-#     elem = create_cg1_quad()
-#     ufl_elem = elem.to_ufl()
-#     print(elem.to_fiat().entity_permutations())
-#     # elem.cell.hasse_diagram(filename="cg1quad.png")
-#     assert (run_test_original(1, "CG", 1, parameters={}, quadrilateral=True) < 1.e-9)
+@pytest.mark.xfail(reason="Issue with quad cell")
+def test_non_tensor_quad():
+    elem = create_cg1_quad()
+    ufl_elem = elem.to_ufl()
+    print(elem.to_fiat().entity_permutations())
+    # elem.cell.hasse_diagram(filename="cg1quad.png")
+    assert (run_test_original(1, "CG", 1, parameters={}, quadrilateral=True) < 1.e-9)
+
+
+def project(U, mesh, func):
+    f = assemble(interpolate(func, U))
+
+    out = Function(U)
+    u = TrialFunction(U)
+    v = TestFunction(U)
+    a = inner(u, v)*dx
+    L = inner(f, v)*dx
+    solve(a == L, out)
+
+    res = sqrt(assemble(dot(out - func, out - func) * dx))
+    return res
 
 
 @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2),