Need to figure out permutation test

ffcx/codegeneration/symbols.py

@@ -138,7 +138,8 @@ def x_component(self, mt):
     def J_component(self, mt):
         """Jacobian component."""
         # FIXME: Add domain number!
-        return L.Symbol(format_mt_name("J", mt), dtype=L.DataType.REAL)
+        return L.Symbol(format_mt_name(
+            f"J{mt.expr.ufl_domain().ufl_id()}", mt), dtype=L.DataType.REAL)
 
     def domain_dof_access(self, dof, component, gdim, num_scalar_dofs, restriction):
         """Domain DOF access."""

ffcx/ir/representationutils.py

@@ -88,7 +88,6 @@ def create_quadrature_points_and_weights(
         if len(edge_types) > 1:
             raise Exception(f"Cell type {cell} not supported for integral type {integral_type}.")
         pts, wts = create_quadrature(edge_types[0].cellname(), degree, rule, elements)
-        print(pts, wts, edge_types)
     elif integral_type in ufl.measure.point_integral_types:
         pts, wts = create_quadrature("vertex", degree, rule, elements)
     elif integral_type == "expression":

test/test_jit_forms.py

@@ -1213,7 +1213,8 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):
     A_ref = tabulate_tensor(ele_type, V_cell_type, V_cell_type, coeffs_ref)
     # Remove the entries for the extra test DOF on the triangle element
     A_ref = A_ref[1:][:]
-
+    print(A_ref)
+    print(A)
     # If the permutation is 1, this means the triangle sees its edge as being flipped
     # relative to the edge's global orientation. Thus the result is the same as swapping
     # cols 1 and 2 and cols 4 and 5 of the reference result.
@@ -1226,14 +1227,15 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):
 
 @pytest.mark.parametrize("dtype", ["float64"])
 @pytest.mark.parametrize("permutation", [[0], [1]])
-def test_mixed_dim_form_codim2(compile_args, dtype, permutation):
+@pytest.mark.parametrize("local_entity_index", [0, 1, 2, 3, 4, 5])
+def test_mixed_dim_form_codim2(compile_args, dtype, permutation, local_entity_index):
     """Test that the local element tensor corresponding to a mixed-dimensional form is correct.
-    The form involves an integral over a facet of the cell. The trial function and a coefficient f
-    are of codim 0. The test function and a coefficient g are of codim 1. We compare against another
+    The form involves an integral over a edge of the cell. The trial function and a coefficient f
+    are of codim 0. The test function and are of codim 2. We compare against another
     form where the test function and g are codim 0 but have the same trace on the facet.
     """
 
-    def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):
+    def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs, l_i):
         "Helper function to create a form and compute the local element tensor"
         V_ele = basix.ufl.element(ele_type, V_cell_type, 2)
         W_ele = basix.ufl.element(ele_type, W_cell_type, 1)
@@ -1251,8 +1253,9 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):
         f = ufl.Coefficient(V)
 
         dl = ufl.Measure("dl", domain=V_domain)
-        forms = [u * f * q * dl]
-        compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms(
+        forms = [u.dx(0) * f * q * dl]
+
+        compiled_forms, module, _ = ffcx.codegeneration.jit.compile_forms(
             forms, options={"scalar_type": dtype}, cffi_extra_compile_args=compile_args
         )
         form0 = compiled_forms[0]
@@ -1262,11 +1265,11 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):
         A = np.zeros((W_ele.dim, V_ele.dim), dtype=dtype)
         w = np.array(coeffs, dtype=dtype)
         c = np.array([], dtype=dtype)
-        facet = np.array([0], dtype=np.intc)
+        edge = np.array([l_i], dtype=np.intc)
         perm = np.array(permutation, dtype=np.uint8)
 
         xdtype = dtype_to_scalar_dtype(dtype)
-        coords = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=xdtype)
+        coords = np.array([[0.0, 0.0, 0.0], [2.5, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0,0.0, 1.3]], dtype=xdtype)
 
         c_type = dtype_to_c_type(dtype)
         c_xtype = dtype_to_c_type(xdtype)
@@ -1277,7 +1280,7 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):
             ffi.cast(f"{c_type}  *", w.ctypes.data),
             ffi.cast(f"{c_type}  *", c.ctypes.data),
             ffi.cast(f"{c_xtype} *", coords.ctypes.data),
-            ffi.cast("int *", facet.ctypes.data),
+            ffi.cast("int *", edge.ctypes.data),
             ffi.cast("uint8_t *", perm.ctypes.data),
         )
 
@@ -1291,29 +1294,55 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):
 
     # Coefficient data
     # f is a quadratic on each edge that is 0 at the vertices and 1 at the midpoint
-    f_data = [0, 0, 0, 1, 1, 1]
-    # g is a linear function along the edge that is 0 at one vertex and 1 at the other
-    g_data = [0, 1]
+    f_data = [5, 6, 7, 1, 2, 3, 4, 5, 6]
+
     # Collect coefficient data
-    coeffs = f_data + g_data
+    coeffs = f_data
 
     # Tabulate the tensor for the mixed-dimensional form
-    A = tabulate_tensor(ele_type, V_cell_type, Vbar_cell_type, coeffs)
+    A = tabulate_tensor(ele_type, V_cell_type, Vbar_cell_type, coeffs, local_entity_index)
 
     # Compare to a reference result. Here, we compare against the same kernel but with
     # the interval element replaced with a triangle.
-    # We create some data for g on the triangle whose trace coincides with g on the interval
-    g_data = [0, 0, 1]
-    coeffs_ref = f_data + g_data
-    A_ref = tabulate_tensor(ele_type, V_cell_type, V_cell_type, coeffs_ref)
-    # Remove the entries for the extra test DOF on the triangle element
-    A_ref = A_ref[1:][:]
+    coeffs_ref = f_data
+    A_ref = tabulate_tensor(ele_type, V_cell_type, V_cell_type, coeffs_ref, local_entity_index)
 
-    # If the permutation is 1, this means the triangle sees its edge as being flipped
-    # relative to the edge's global orientation. Thus the result is the same as swapping
-    # cols 1 and 2 and cols 4 and 5 of the reference result.
+    # Find the local indices of the dofs that are on the edge
+    local_index_to_slice = {0: [2,3],
+                            1: [1,3],
+                            2: [1,2],
+                            3: [0,3],
+                            4: [0,2],
+                            5: [0,1]}
+
+
+    A_ref = A_ref[local_index_to_slice[local_entity_index]]
+
+    print(A_ref.T)
+    edge_flips = {0: [[2,3]],#,[5,6],[7,8]],
+                      1: [[1,3],[4,6],[7,9]],
+                      2: [[1,2],[4,5],[8,9]],
+                        3: [[0,3],[5,9],[4,8]],
+                        4: [[0,2],[4,7],[6,9]],
+                        5: [[0,1],[5,7],[6,8]]
+                }
+
+
+    # We only check a single edge flip (as it is 
     if permutation[0] == 1:
-        A_ref[:, [1, 2]] = A_ref[:, [2, 1]]
-        A_ref[:, [4, 5]] = A_ref[:, [5, 4]]
+        if local_entity_index == 0:
+            for flips in edge_flips[local_entity_index]:
+                A_ref[:,flips] = A_ref[:,flips[::-1]]
+        A_ref[:] = A_ref[::-1]    
+    # # If the permutation is 1, this means the triangle sees its edge as being flipped
+    # # relative to the edge's global orientation. Thus the result is the same as swapping
+    # # cols 1 and 2 and cols 4 and 5 of the reference result.
+    # if permutation[0] == 1:
+    #     A_ref[:, [1, 2]] = A_ref[:, [2, 1]]
+    #     A_ref[:, [4, 5]] = A_ref[:, [5, 4]]
+    # print(A.T)
+    print()
+    print(A_ref.T)
+    print(A.T)
 
     assert np.allclose(A, A_ref)