Indiamai/addition

Makefile

@@ -26,6 +26,12 @@ tests:
 	@echo "    Running all tests"
 	@FIREDRAKE_USE_FUSE=1 python3 -m coverage run -p -m pytest -rx test
 
+mini_tests:
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_2d_examples_docs.py
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_1d
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_orientations.py::test_surface_vec_rt
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_projection_convergence_3d\[construct_tet_ned-N1curl-1-0.8\]
+
 coverage:
 	@python3 -m coverage combine
 	@python3 -m coverage report -m 

fuse/dof.py

@@ -347,14 +347,13 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non
             self.pairing = self.pairing.add_entity(cell)
         if self.target_space is None:
             self.target_space = space
-        if self.id is None and overall_id is not None:
+        if overall_id is not None:
             self.id = overall_id
-        if self.sub_id is None and generator_id is not None:
+        if generator_id is not None:
             self.sub_id = generator_id
 
     def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()):
         # TODO deriv dict needs implementing (currently {})
-        print(value_shape)
         return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree, value_shape), {}, str(self))
 
     def to_quadrature(self, arg_degree, value_shape):

fuse/spaces/element_sobolev_spaces.py

@@ -1,6 +1,7 @@
 from functools import total_ordering
 from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2
 
+
 @total_ordering
 class ElementSobolevSpace(object):
     """
@@ -54,7 +55,7 @@ def __init__(self, cell):
 
     def __repr__(self):
         return "H1"
-    
+
     def to_ufl(self):
         return H1
 

fuse/spaces/polynomial_spaces.py

@@ -96,9 +96,7 @@ def __mul__(self, x):
         the sympy object on the right. This is due to Sympy's implementation of __mul__ not
         passing to this handler as it should.
         """
-        if isinstance(x, sp.Symbol):
-            return ConstructedPolynomialSpace([x], [self])
-        elif isinstance(x, sp.Matrix):
+        if isinstance(x, sp.Symbol) or isinstance(x, sp.Expr) or isinstance(x, sp.Matrix):
             return ConstructedPolynomialSpace([x], [self])
         else:
             raise TypeError(f'Cannot multiply a PolySpace with {type(x)}')
@@ -177,8 +175,6 @@ def to_ON_polynomial_set(self, ref_el):
 
         # otherwise have to work on this through tabulation
 
-        Q = create_quadrature(ref_el, 2 * (k + 1))
-        Qpts, Qwts = Q.get_points(), Q.get_weights()
         weighted_sets = []
 
         for (s, w) in zip(self.spaces, self.weights):
@@ -196,6 +192,8 @@ def to_ON_polynomial_set(self, ref_el):
                 else:
                     vec = True
                 w_deg = max_deg_sp_expr(w)
+                Q = create_quadrature(ref_el, 2 * (k + w_deg + 1))
+                Qpts, Qwts = Q.get_points(), Q.get_weights()
                 Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, shape, scale="orthonormal")
                 # vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal")
 
@@ -206,11 +204,11 @@ def to_ON_polynomial_set(self, ref_el):
                 if s.set_shape or vec:
                     scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :]
                 else:
-                    # breakpoint()
                     scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :]
                     scaled_at_Qpts = scaled_at_Qpts.squeeze()
                 PkHw_coeffs = np.dot(np.multiply(scaled_at_Qpts, Qwts), Pkpw_at_Qpts.T)
-                # breakpoint()
+                if len(PkHw_coeffs.shape) == 1:
+                    PkHw_coeffs = PkHw_coeffs.reshape(1, -1)
                 weighted_sets.append(polynomial_set.PolynomialSet(ref_el,
                                                                   space.degree + w_deg,
                                                                   space.degree + w_deg,

fuse/triples.py

@@ -421,7 +421,7 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set):
                         ent_dofs = entity_associations[dim][e_id][dof_gen]
                         ent_dofs_ids = np.array([self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs], dtype=int)
                         # (dof_gen, ent_dofs)
-                        total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if ed.id not in total_ent_dof_ids]
+                        total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if self.dof_id_to_fiat_id[ed.id] not in total_ent_dof_ids]
                         # dof_idx = [total_ent_dof_ids.index(id) for id in ent_dofs_ids]
                         dof_gen_class = ent_dofs[0].generation
 
@@ -537,6 +537,21 @@ def reverse_dof_perms(self, matrices):
                 reversed_mats[dim][e_id] = perms_copy
         return reversed_mats
 
+    def __add__(self, other):
+        """ Construct a new element triple by combining the degrees of freedom
+        This implementation does not make assertions about the properties
+        of the resulting element.
+
+        Elements being adding must be defined over the same cell and have the same
+        value shape and mapping"""
+        assert self.cell == other.cell
+        assert self.spaces[0].set_shape == other.spaces[0].set_shape
+        assert str(self.spaces[1]) == str(other.spaces[1])
+
+        spaces = (self.spaces[0] + other.spaces[0], self.spaces[1], max([self.spaces[2], other.spaces[2]]))
+
+        return ElementTriple(self.cell, spaces, self.DOFGenerator + other.DOFGenerator)
+
     def _to_dict(self):
         o_dict = {"cell": self.cell, "spaces": self.spaces, "dofs": self.DOFGenerator}
         return o_dict
@@ -574,21 +589,20 @@ def num_dofs(self):
         return self.dof_numbers
 
     def generate(self, cell, space, id_counter):
-        if self.ls is None:
-            self.ls = []
-            for l_g in self.x:
-                i = 0
-                for g in self.g1.members():
-                    generated = l_g(g)
-                    if not isinstance(generated, list):
-                        generated = [generated]
-                    for dof in generated:
-                        dof.add_context(self, cell, space, g, id_counter, i)
-                        id_counter += 1
-                        i += 1
-                    self.ls.extend(generated)
-            self.dof_numbers = len(self.ls)
-            self.dof_ids = [dof.id for dof in self.ls]
+        self.ls = []
+        for l_g in self.x:
+            i = 0
+            for g in self.g1.members():
+                generated = l_g(g)
+                if not isinstance(generated, list):
+                    generated = [generated]
+                for dof in generated:
+                    dof.add_context(self, cell, space, g, id_counter, i)
+                    id_counter += 1
+                    i += 1
+                self.ls.extend(generated)
+        self.dof_numbers = len(self.ls)
+        self.dof_ids = [dof.id for dof in self.ls]
         return self.ls
 
     def make_entity_ids(self):

fuse/utils.py

@@ -67,11 +67,10 @@ def max_deg_sp_expr(sp_expr):
         for comp in sp_expr:
             # only compute degree if component is a polynomial
             if sp.sympify(comp).as_poly():
-                degs += [sp.sympify(comp).as_poly().degree()]
+                degs += [sp.sympify(comp).as_poly().total_degree()]
     else:
         if sp.sympify(sp_expr).as_poly():
-            degs += [sp.sympify(sp_expr).as_poly().degree()]
-
+            degs += [sp.sympify(sp_expr).as_poly().total_degree()]
     return max(degs)
 
 

test/test_3d_examples_docs.py

@@ -327,7 +327,7 @@ def test_tet_rt2():
     ls = rt2.generate()
     # TODO make this a proper test
     for dof in ls:
-        print(dof.to_quadrature(1))
+        print(dof.to_quadrature(1, value_shape=(2,)))
     rt2.to_fiat()
 
 

test/test_algebra.py

@@ -0,0 +1,36 @@
+from fuse import *
+from firedrake import *
+import numpy as np
+import sympy as sp
+from test_convert_to_fiat import create_cg2_tri, construct_cg3
+
+
+def construct_bubble(cell=None):
+    if cell is None:
+        cell = polygon(3)
+    x = sp.Symbol("x")
+    y = sp.Symbol("y")
+    f = (3*np.sqrt(3)/4)*(y + np.sqrt(3)/3)*(np.sqrt(3)*x + y - 2*np.sqrt(3)/3)*(-np.sqrt(3)*x + y - 2*np.sqrt(3)/3)
+    space = PolynomialSpace(3).restrict(0, 0)*f
+    xs = [DOF(DeltaPairing(), PointKernel((0, 0)))]
+    bubble = ElementTriple(cell, (space, CellL2, L2), DOFGenerator(xs, S1, S1))
+    return bubble
+
+
+def test_bubble():
+    mesh = UnitTriangleMesh()
+    x = SpatialCoordinate(mesh)
+
+    tri = polygon(3)
+    bub = construct_bubble(tri)
+    cg2 = create_cg2_tri(tri)
+    p2b3 = bub + cg2
+    V = FunctionSpace(mesh, p2b3.to_ufl())
+    W = FunctionSpace(mesh, construct_cg3().to_ufl())
+
+    bubble_func = 27*x[0]*x[1]*(1-x[0]-x[1])
+    u = project(bubble_func, V)
+    exact = Function(W)
+    exact.interpolate(bubble_func, W)
+    # make sure that these are the same
+    assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14

test/test_convert_to_fiat.py

@@ -358,7 +358,7 @@ def test_entity_perms(elem_gen, cell):
 
 @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg1, "CG", 1),
                                                     (create_dg1, "DG", 1),
-                                                    (construct_dg0_integral, "DG", 0),
+                                                    pytest.param(construct_dg0_integral, "DG", 0, marks=pytest.mark.xfail(reason='Passes locally, fails in CI, probably same as dg2')),
                                                     (construct_dg1_integral, "DG", 1),
                                                     (construct_dg2_integral, "DG", 2),
                                                     pytest.param(create_dg2, "DG", 2, marks=pytest.mark.xfail(reason='Need to update TSFC in CI')),

test/test_dofs.py

@@ -220,6 +220,6 @@ def test_generate_quadrature():
         print("fiat", d.pt_dict)
     print()
     for d in elem.generate():
-        print("fuse", d.to_quadrature(degree))
+        print("fuse", d.to_quadrature(degree, value_shape=(2,)))
 
     elem.to_fiat()

test/test_perms.py

@@ -1,9 +1,8 @@
 from fuse import *
-from test_convert_to_fiat import create_cg1, create_dg1, create_cg2
+from test_convert_to_fiat import create_cg1, create_dg1, create_cg2, construct_nd
 from test_2d_examples_docs import construct_cg3
 import pytest
 import numpy as np
-import sympy as sp
 
 vert = Point(0)
 edge = Point(1, [Point(0), Point(0)], vertex_num=2)
@@ -48,24 +47,7 @@ def test_basic_perms(cell):
 
 @pytest.mark.parametrize("cell", [tri])
 def test_nd_perms(cell):
-    deg = 1
-    edge = cell.edges(get_class=True)[0]
-    x = sp.Symbol("x")
-    y = sp.Symbol("y")
-
-    xs = [DOF(L2Pairing(), PolynomialKernel(edge.basis_vectors()[0]))]
-    dofs = DOFGenerator(xs, S1, S2)
-    int_ned = ElementTriple(edge, (P1, CellHCurl, C0), dofs)
-
-    xs = [immerse(cell, int_ned, TrHCurl)]
-    tri_dofs = DOFGenerator(xs, C3, S3)
-
-    M = sp.Matrix([[y, -x]])
-    vec_Pk = PolynomialSpace(deg - 1, set_shape=True)
-    Pk = PolynomialSpace(deg - 1)
-    nd = vec_Pk + (Pk.restrict(deg - 2, deg - 1))*M
-
-    ned = ElementTriple(cell, (nd, CellHCurl, C0), [tri_dofs])
+    ned = construct_nd(cell)
     ned.to_fiat()
     for i, mat in ned.matrices[2][0].items():
         print(i)

test/test_tensor_prod.py

@@ -36,10 +36,10 @@ def mass_solve(U):
 
 
 @pytest.mark.parametrize("generator1, generator2, code1, code2, deg1, deg2",
-                          [(construct_cg1, construct_cg1, "CG", "CG", 1, 1),
-                           (construct_dg1, construct_dg1, "DG", "DG", 1, 1),
-                           (construct_dg1, construct_cg1, "DG", "CG", 1, 1),
-                           (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)])
+                         [(construct_cg1, construct_cg1, "CG", "CG", 1, 1),
+                          (construct_dg1, construct_dg1, "DG", "DG", 1, 1),
+                          (construct_dg1, construct_cg1, "DG", "CG", 1, 1),
+                          (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)])
 def test_ext_mesh(generator1, generator2, code1, code2, deg1, deg2):
     m = UnitIntervalMesh(2)
     mesh = ExtrudedMesh(m, 2)
@@ -137,3 +137,59 @@ def test_flattening(A, B, res):
     else:
         cell = tensor_cell.flatten()
         cell.construct_fuse_rep()
+
+
+# def test_trace_galerkin_projection():
+#     mesh = UnitSquareMesh(10, 10, quadrilateral=True)
+
+#     x, y = SpatialCoordinate(mesh)
+#     A = construct_cg1()
+#     B = construct_dg0_integral()
+#     elem = tensor_product(A, B)
+#     elem = elem.flatten()
+
+#     # Define the Trace Space
+#     T = FunctionSpace(mesh, elem.to_ufl())
+
+#     # Define trial and test functions
+#     lambdar = TrialFunction(T)
+#     gammar = TestFunction(T)
+
+#     # Define right hand side function
+
+#     V = FunctionSpace(mesh, "CG", 1)
+#     f = Function(V)
+#     f.interpolate(cos(x*pi*2)*cos(y*pi*2))
+
+#     # Construct bilinear form
+#     a = inner(lambdar, gammar) * ds + inner(lambdar('+'), gammar('+')) * dS
+
+#     # Construct linear form
+#     l = inner(f, gammar) * ds + inner(f('+'), gammar('+')) * dS
+
+#     # Compute the solution
+#     t = Function(T)
+#     solve(a == l, t, solver_parameters={'ksp_rtol': 1e-14})
+
+#     # Compute error in trace norm
+#     trace_error = sqrt(assemble(FacetArea(mesh)*inner((t - f)('+'), (t - f)('+')) * dS))
+
+#     assert trace_error < 1e-13
+
+# def test_hdiv():
+#     np.set_printoptions(linewidth=90, precision=4, suppress=True)
+#     m = UnitIntervalMesh(2)
+#     mesh = ExtrudedMesh(m, 2)
+#     CG_1 = FiniteElement("CG", interval, 1)
+#     DG_0 = FiniteElement("DG", interval, 0)
+#     P1P0 = TensorProductElement(CG_1, DG_0)
+#     RT_horiz = HDivElement(P1P0)
+#     P0P1 = TensorProductElement(DG_0, CG_1)
+#     RT_vert = HDivElement(P0P1)
+#     elt = RT_horiz + RT_vert
+#     V = FunctionSpace(mesh, elt)
+#     tabulation = V.finat_element.fiat_equivalent.tabulate(0, [(0, 0), (1, 0)])
+#     for ent, arr in tabulation.items():
+#         print(ent)
+#         for comp in arr:
+#             print(comp[0], comp[1])