From dcf329b4e4256a98839896bbff058a0fb4207b49 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 16 Dec 2025 15:05:45 +0100 Subject: [PATCH 1/8] Add test of quadrature element and check value shape --- tests/test_interpolation.py | 33 +++++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index f4b2cb5..d13d04b 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -5,6 +5,7 @@ import pytest import ufl import numpy as np +import basix.ufl @pytest.mark.skipif( @@ -20,8 +21,10 @@ ], ) @pytest.mark.parametrize("use_petsc", [True, False]) -@pytest.mark.parametrize("degree", [1, 3, 5]) -def test_interpolation_matrix(use_petsc, cell_type, degree): +@pytest.mark.parametrize("degree", [2, 4]) +@pytest.mark.parametrize("out_family", ["Lagrange", "DG", "Quadrature"]) +@pytest.mark.parametrize("value_shape", [(), (2,), (2, 3)]) +def test_interpolation_matrix(use_petsc, cell_type, degree, out_family, value_shape): if use_petsc: pytest.importorskip("petsc4py") @@ -33,14 +36,24 @@ def test_interpolation_matrix(use_petsc, cell_type, degree): else: raise ValueError("Unsupported cell type") - V = dolfinx.fem.functionspace(mesh, ("DG", degree)) - Q = dolfinx.fem.functionspace(mesh, ("Lagrange", degree)) + V = dolfinx.fem.functionspace(mesh, ("DG", degree, value_shape)) + if out_family == "Quadrature": + el = basix.ufl.quadrature_element(mesh.basix_cell(), degree=degree, value_shape=value_shape) + else: + el = (out_family, degree, value_shape) + Q = dolfinx.fem.functionspace(mesh, el) + + def f(x): + scalar_val = x[0] ** degree + x[1] if tdim == 2 else x[0] + x[1] + x[2] ** degree + vs = int(np.prod(value_shape)) + return np.tile(scalar_val, vs).reshape(vs, -1) u = dolfinx.fem.Function(V) - u.interpolate(lambda x: x[0] ** degree + x[1] if tdim == 2 else x[0] + x[1] + x[2] ** degree) + u.interpolate(f) q = dolfinx.fem.Function(Q) expr = ufl.TrialFunction(V) + if use_petsc: A = scifem.interpolation.petsc_interpolation_matrix(expr, Q) A.mult(u.x.petsc_vec, q.x.petsc_vec) @@ -59,7 +72,15 @@ def test_interpolation_matrix(use_petsc, cell_type, degree): q.x.scatter_forward() q_ref = dolfinx.fem.Function(Q) - q_ref.interpolate(u) + if out_family == "Quadrature": + try: + ip = Q.element.interpolation_points() + except TypeError: + ip = Q.element.interpolation_points + u_expr = dolfinx.fem.Expression(u, ip) + q_ref.interpolate(u_expr) + else: + q_ref.interpolate(u) np.testing.assert_allclose(q.x.array, q_ref.x.array, rtol=1e-12, atol=1e-13) From 80b1da79f563dff5a11ebd89fb2337c809ad9af0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 16 Dec 2025 15:17:18 +0100 Subject: [PATCH 2/8] Add fixes for quadrature elements --- src/scifem/interpolation.py | 35 ++++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index ed32a33..91031f0 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -55,8 +55,6 @@ def prepare_interpolation_data( array_evaluated = compiled_expr.eval(mesh, np.arange(num_cells, dtype=np.int32)) assert np.prod(Q.value_shape) == np.prod(expr.ufl_shape) - im = Q.element.basix_element.interpolation_matrix - # Get data as (num_cells*num_points,1, expr_shape, num_test_basis_functions*test_block_size) expr_size = int(np.prod(expr.ufl_shape)) array_evaluated = array_evaluated.reshape( @@ -73,19 +71,29 @@ def prepare_interpolation_data( num_cells * num_points, mesh.geometry.dim, mesh.topology.dim ) - Q_vs = Q.element.basix_element.value_size + try: + Q_vs = Q.element.basix_element.value_size + except RuntimeError: + Q_vs = 1 # If we do not have a basix element, assume value size is 1 + new_array = np.zeros( (num_cells * num_points, Q.dofmap.bs * Q_vs, V.dofmap.bs * V.dofmap.dof_layout.num_dofs), dtype=np.float64, ) - for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): - for q in range(Q.dofmap.bs): - new_array[:, q * Q_vs : (q + 1) * Q_vs, i] = Q.element.basix_element.pull_back( - array_evaluated[:, :, q * Q_vs : (q + 1) * Q_vs, i], jacs, detJs, Ks - ).reshape(num_cells * num_points, Q_vs) - new_array = new_array.reshape( - num_cells, num_points, Q.dofmap.bs * Q_vs, V.dofmap.bs * V.dofmap.dof_layout.num_dofs - ) + if not isinstance(Q.ufl_element().pullback, ufl.pullback.IdentityPullback): + for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): + for q in range(Q.dofmap.bs): + new_array[:, q * Q_vs : (q + 1) * Q_vs, i] = Q.element.basix_element.pull_back( + array_evaluated[:, :, q * Q_vs : (q + 1) * Q_vs, i], jacs, detJs, Ks + ).reshape(num_cells * num_points, Q_vs) + new_array = new_array.reshape( + num_cells, num_points, Q.dofmap.bs * Q_vs, V.dofmap.bs * V.dofmap.dof_layout.num_dofs + ) + else: + new_array = array_evaluated.reshape( + num_cells, num_points, Q.dofmap.bs * Q_vs, V.dofmap.bs * V.dofmap.dof_layout.num_dofs + ) + interpolated_matrix = np.zeros( ( num_cells, @@ -95,6 +103,11 @@ def prepare_interpolation_data( dtype=np.float64, ) + if Q.element.interpolation_ident: + im = np.eye(Q.element.interpolation_points.shape[0]) + else: + im = Q.element.basix_element.interpolation_matrix + for c in range(num_cells): for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): tmp_array = np.zeros((int(num_points), Q.dofmap.bs * Q_vs), dtype=np.float64) From 259d4e7a3c48387162dde4ee709d2ea81b05f030 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 16 Dec 2025 16:57:10 +0000 Subject: [PATCH 3/8] Add vectorized version for nice spaces --- src/scifem/interpolation.py | 36 ++++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index 91031f0..020676b 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -104,22 +104,30 @@ def prepare_interpolation_data( ) if Q.element.interpolation_ident: - im = np.eye(Q.element.interpolation_points.shape[0]) + # Smart vectorized version with identity mapping + if Q.dofmap.bs == 1: + interpolated_matrix = new_array.transpose(0, 2, 1, 3).reshape(new_array.shape[0], new_array.shape[1] * new_array.shape[2], new_array.shape[3]) + else: + i_scalar = new_array.transpose(0, 2, 1, 3) + interpolated_matrix = np.zeros((new_array.shape[0], new_array.shape[1] * new_array.shape[2], new_array.shape[3])) + for q in range(Q.dofmap.bs): + interpolated_matrix[:, q::Q.dofmap.bs, :] = i_scalar[:, q, :, :] + else: + # Tedious non-identity version im = Q.element.basix_element.interpolation_matrix - - for c in range(num_cells): - for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): - tmp_array = np.zeros((int(num_points), Q.dofmap.bs * Q_vs), dtype=np.float64) - for p in range(num_points): - tmp_array[p] = new_array[c, p, :, i] - if Q.dofmap.bs == 1: - interpolated_matrix[c, :, i] = (im @ tmp_array.T.flatten()).flatten() - else: - for q in range(Q.dofmap.bs): - interpolated_matrix[c, q :: Q.dofmap.bs, i] = ( - im @ tmp_array.T[q].flatten() - ).flatten() + for c in range(num_cells): + for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): + tmp_array = np.zeros((int(num_points), Q.dofmap.bs * Q_vs), dtype=np.float64) + for p in range(num_points): + tmp_array[p] = new_array[c, p, :, i] + if Q.dofmap.bs == 1: + interpolated_matrix[c, :, i] = (im @ tmp_array.T.flatten()).flatten() + else: + for q in range(Q.dofmap.bs): + interpolated_matrix[c, q :: Q.dofmap.bs, i] = ( + im @ tmp_array.T[q].flatten() + ).flatten() # Apply dof transformation to each column (using Piopla maps) mesh.topology.create_entity_permutations() if Q.element.needs_dof_transformations: From 89268e429330e74cf6824a40cb7e3c7d0786d1a1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 16 Dec 2025 16:57:26 +0000 Subject: [PATCH 4/8] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/scifem/interpolation.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index 020676b..4ce0a66 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -106,12 +106,16 @@ def prepare_interpolation_data( if Q.element.interpolation_ident: # Smart vectorized version with identity mapping if Q.dofmap.bs == 1: - interpolated_matrix = new_array.transpose(0, 2, 1, 3).reshape(new_array.shape[0], new_array.shape[1] * new_array.shape[2], new_array.shape[3]) + interpolated_matrix = new_array.transpose(0, 2, 1, 3).reshape( + new_array.shape[0], new_array.shape[1] * new_array.shape[2], new_array.shape[3] + ) else: i_scalar = new_array.transpose(0, 2, 1, 3) - interpolated_matrix = np.zeros((new_array.shape[0], new_array.shape[1] * new_array.shape[2], new_array.shape[3])) + interpolated_matrix = np.zeros( + (new_array.shape[0], new_array.shape[1] * new_array.shape[2], new_array.shape[3]) + ) for q in range(Q.dofmap.bs): - interpolated_matrix[:, q::Q.dofmap.bs, :] = i_scalar[:, q, :, :] + interpolated_matrix[:, q :: Q.dofmap.bs, :] = i_scalar[:, q, :, :] else: # Tedious non-identity version From 18668427900810283e3681a83c21dfcd67075baa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 16 Dec 2025 17:07:12 +0000 Subject: [PATCH 5/8] Final optimization, compute jacobians only in case of pullback required --- src/scifem/interpolation.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index 4ce0a66..e4d8728 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -60,16 +60,6 @@ def prepare_interpolation_data( array_evaluated = array_evaluated.reshape( num_cells * q_points.shape[0], 1, expr_size, V.dofmap.bs * V.dofmap.dof_layout.num_dofs ) - jacobian = dolfinx.fem.Expression(ufl.Jacobian(mesh), q_points) - detJ = dolfinx.fem.Expression(ufl.JacobianDeterminant(mesh), q_points) - K = dolfinx.fem.Expression(ufl.JacobianInverse(mesh), q_points) - jacs = jacobian.eval(mesh, np.arange(num_cells, dtype=np.int32)).reshape( - num_cells * num_points, mesh.geometry.dim, mesh.topology.dim - ) - detJs = detJ.eval(mesh, np.arange(num_cells, dtype=np.int32)).flatten() - Ks = K.eval(mesh, np.arange(num_cells, dtype=np.int32)).reshape( - num_cells * num_points, mesh.geometry.dim, mesh.topology.dim - ) try: Q_vs = Q.element.basix_element.value_size @@ -81,6 +71,17 @@ def prepare_interpolation_data( dtype=np.float64, ) if not isinstance(Q.ufl_element().pullback, ufl.pullback.IdentityPullback): + jacobian = dolfinx.fem.Expression(ufl.Jacobian(mesh), q_points) + detJ = dolfinx.fem.Expression(ufl.JacobianDeterminant(mesh), q_points) + K = dolfinx.fem.Expression(ufl.JacobianInverse(mesh), q_points) + jacs = jacobian.eval(mesh, np.arange(num_cells, dtype=np.int32)).reshape( + num_cells * num_points, mesh.geometry.dim, mesh.topology.dim + ) + detJs = detJ.eval(mesh, np.arange(num_cells, dtype=np.int32)).flatten() + Ks = K.eval(mesh, np.arange(num_cells, dtype=np.int32)).reshape( + num_cells * num_points, mesh.geometry.dim, mesh.topology.dim + ) + for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): for q in range(Q.dofmap.bs): new_array[:, q * Q_vs : (q + 1) * Q_vs, i] = Q.element.basix_element.pull_back( From ae1a2a87e78fbf6c8daa21f77e2dd7acec54718a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 16 Dec 2025 17:12:13 +0000 Subject: [PATCH 6/8] Avoid repeat structures in testing function for value size --- tests/test_interpolation.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index d13d04b..71d23f1 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -46,7 +46,10 @@ def test_interpolation_matrix(use_petsc, cell_type, degree, out_family, value_sh def f(x): scalar_val = x[0] ** degree + x[1] if tdim == 2 else x[0] + x[1] + x[2] ** degree vs = int(np.prod(value_shape)) - return np.tile(scalar_val, vs).reshape(vs, -1) + f_rep = np.tile(scalar_val, vs).reshape(vs, -1) + for i in range(vs): + f_rep[i] += np.pi * (i+1) + return f_rep u = dolfinx.fem.Function(V) u.interpolate(f) From eb7c8d3ad42ffb51c8ddb6e73c02bbfd88950fc1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 16 Dec 2025 17:12:26 +0000 Subject: [PATCH 7/8] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_interpolation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 71d23f1..3d48732 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -48,7 +48,7 @@ def f(x): vs = int(np.prod(value_shape)) f_rep = np.tile(scalar_val, vs).reshape(vs, -1) for i in range(vs): - f_rep[i] += np.pi * (i+1) + f_rep[i] += np.pi * (i + 1) return f_rep u = dolfinx.fem.Function(V) From df11038d06fbc6334376c372b0e03961f95f2af7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 17 Dec 2025 10:05:34 +0100 Subject: [PATCH 8/8] Clarify usage and error handling of q-els --- src/scifem/interpolation.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index e4d8728..effb244 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -61,15 +61,26 @@ def prepare_interpolation_data( num_cells * q_points.shape[0], 1, expr_size, V.dofmap.bs * V.dofmap.dof_layout.num_dofs ) + # Check if we are dealing with a quadrature element or not. + # They do not have a complete DOLFINx API, which makes them tricky to use. try: - Q_vs = Q.element.basix_element.value_size + basix_el = Q.element.basix_element + Q_vs = basix_el.value_size + pull_back = basix_el.pull_back + im = basix_el.interpolation_matrix except RuntimeError: Q_vs = 1 # If we do not have a basix element, assume value size is 1 + assert isinstance(Q.ufl_element().pullback, ufl.pullback.IdentityPullback) + pull_back = lambda x: None + assert Q.element.interpolation_ident + im = None new_array = np.zeros( (num_cells * num_points, Q.dofmap.bs * Q_vs, V.dofmap.bs * V.dofmap.dof_layout.num_dofs), dtype=np.float64, ) + + # Check if pullback is identity, then we can skip this step if not isinstance(Q.ufl_element().pullback, ufl.pullback.IdentityPullback): jacobian = dolfinx.fem.Expression(ufl.Jacobian(mesh), q_points) detJ = dolfinx.fem.Expression(ufl.JacobianDeterminant(mesh), q_points) @@ -84,7 +95,7 @@ def prepare_interpolation_data( for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): for q in range(Q.dofmap.bs): - new_array[:, q * Q_vs : (q + 1) * Q_vs, i] = Q.element.basix_element.pull_back( + new_array[:, q * Q_vs : (q + 1) * Q_vs, i] = pull_back( array_evaluated[:, :, q * Q_vs : (q + 1) * Q_vs, i], jacs, detJs, Ks ).reshape(num_cells * num_points, Q_vs) new_array = new_array.reshape( @@ -103,7 +114,8 @@ def prepare_interpolation_data( ), dtype=np.float64, ) - + # Check if interpolation matrix of dual operator is identity, then we can use a vectorized + # version of this step if Q.element.interpolation_ident: # Smart vectorized version with identity mapping if Q.dofmap.bs == 1: @@ -120,7 +132,6 @@ def prepare_interpolation_data( else: # Tedious non-identity version - im = Q.element.basix_element.interpolation_matrix for c in range(num_cells): for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): tmp_array = np.zeros((int(num_points), Q.dofmap.bs * Q_vs), dtype=np.float64)