From 5fa9d305bf862921d2fb0e821a82998740da4f51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 5 Dec 2024 17:45:14 +0000 Subject: [PATCH] Need to figure out permutation test --- ffcx/codegeneration/symbols.py | 3 +- ffcx/ir/representationutils.py | 1 - test/test_jit_forms.py | 81 +++++++++++++++++++++++----------- 3 files changed, 57 insertions(+), 28 deletions(-) diff --git a/ffcx/codegeneration/symbols.py b/ffcx/codegeneration/symbols.py index 236e5d983..5fe006337 100644 --- a/ffcx/codegeneration/symbols.py +++ b/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.""" diff --git a/ffcx/ir/representationutils.py b/ffcx/ir/representationutils.py index a1725367a..68e705038 100644 --- a/ffcx/ir/representationutils.py +++ b/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": diff --git a/test/test_jit_forms.py b/test/test_jit_forms.py index 6187027cf..1df7e2d49 100644 --- a/test/test_jit_forms.py +++ b/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)