diff --git a/CITATION.cff b/CITATION.cff
index c64d066..598c487 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -1,6 +1,6 @@
 abstract: "<p>Major improvements for interpolation matrix.</p>
 
-  <p><strong>Full Changelog</strong>: https://github.com/scientificcomputing/scifem/compare/v0.16.0...v0.15.2</p>"
+  <p><strong>Full Changelog</strong>: https://github.com/scientificcomputing/scifem/compare/v0.16.1...v0.16.0</p>"
 authors:
   - affiliation: Simula Research Laboratoy
     family-names: Henrik Finsberg
@@ -10,7 +10,7 @@ cff-version: 1.2.0
 date-released: "2025-12-17"
 doi: 10.5281/zenodo.13784777
 license: MIT
-repository-code: https://github.com/scientificcomputing/scifem/tree/v0.16.0
-title: "scientificcomputing/scifem: v0.16.0"
+repository-code: https://github.com/scientificcomputing/scifem/tree/v0.16.1
+title: "scientificcomputing/scifem: v0.16.1"
 type: software
-version: v0.16.0
+version: v0.16.1
diff --git a/pyproject.toml b/pyproject.toml
index eb157fb..edce138 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,7 +10,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "scifem"
-version = "0.16.0"
+version = "0.16.1"
 description = "Scientific tools for finite element methods"
 readme = "README.md"
 requires-python = ">=3.10"
@@ -138,7 +138,7 @@ tag = true
 sign_tags = false
 tag_name = "v{new_version}"
 tag_message = "Bump version: {current_version} → {new_version}"
-current_version = "0.16.0"
+current_version = "0.16.1"
 
 
 [[tool.bumpversion.files]]
diff --git a/src/scifem/__init__.py b/src/scifem/__init__.py
index 6f2cd3a..b1a0a24 100644
--- a/src/scifem/__init__.py
+++ b/src/scifem/__init__.py
@@ -1,4 +1,5 @@
 from __future__ import annotations
+from importlib.metadata import metadata
 
 import dolfinx
 import numpy as np
@@ -22,8 +23,20 @@
 from .eval import evaluate_function
 from .interpolation import interpolation_matrix, prepare_interpolation_data
 
+meta = metadata("scifem")
+__version__ = meta["Version"]
+__author__ = meta.get("Author", "")
+__license__ = meta.get("License", "MIT")
+__email__ = meta["Author-email"]
+__program_name__ = meta["Name"]
 
 __all__ = [
+    "meta",
+    "__version__",
+    "__author__",
+    "__license__",
+    "__email__",
+    "__program_name__",
     "PointSource",
     "assemble_scalar",
     "create_space_of_simple_functions",
diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py
index effb244..89aad00 100644
--- a/src/scifem/interpolation.py
+++ b/src/scifem/interpolation.py
@@ -14,7 +14,9 @@
 
 
 def prepare_interpolation_data(
-    expr: ufl.core.expr.Expr, Q: dolfinx.fem.FunctionSpace
+    expr: ufl.core.expr.Expr,
+    Q: dolfinx.fem.FunctionSpace,
+    interpolation_entities: npt.NDArray[np.int32] | None = None,
 ) -> npt.NDArray[np.inexact]:
     """Convenience function for preparing data required for assembling the interpolation matrix
 
@@ -30,29 +32,60 @@ def prepare_interpolation_data(
     Args:
         expr: The UFL expression containing a trial function from space `V`
         Q: Output interpolation space
+        interpolation_entities: Entities of the domain of the input space `V` that one
+            should evaluate the `expr` at. If not provided, it is assumed that
+            we are integrating over all cells in `V` and that `Q` is defined on the same grid.
     Returns:
         Interpolation data per cell, as an numpy array.
     """
     if np.issubdtype(dolfinx.default_scalar_type, np.complexfloating):
         raise NotImplementedError("No complex support")
 
+    # Extract argument from expr (in V)
+    arguments = ufl.algorithms.extract_arguments(expr)
+    assert len(arguments) == 1
+    V = arguments[0].ufl_function_space()
+
+    mesh = V.mesh
+
+    if Q.mesh.topology.dim == V.mesh.topology.dim:
+        if interpolation_entities is None:
+            tdim = mesh.topology.dim
+            num_cells = mesh.topology.index_map(tdim).size_local
+            interpolation_entities = np.arange(num_cells, dtype=np.int32)
+        else:
+            if (ndim := interpolation_entities.ndim) != 1:
+                raise ValueError(
+                    f"Interpolation entities has wrong input shape, should be 1D, got {ndim}"
+                )
+            num_cells = len(interpolation_entities)
+    elif Q.mesh.topology.dim == V.mesh.topology.dim - 1:
+        if interpolation_entities is None:
+            raise ValueError(
+                "For integration onto a submesh of codim 1,"
+                + "the integration entities has to be provided"
+            )
+        else:
+            if (ndim := interpolation_entities.ndim) != 2:
+                raise ValueError(
+                    f"Interpolation entities has wrong input shape, should be 2D, got {ndim}"
+                )
+            num_cells = interpolation_entities.shape[0]
+    else:
+        raise RuntimeError("Only codim-1 interpolation matrices can be defined")
+
+    # Extract quadrature points for expression (in Q space)
     try:
         q_points = Q.element.interpolation_points()
     except TypeError:
         q_points = Q.element.interpolation_points
 
-    arguments = ufl.algorithms.extract_arguments(expr)
-    assert len(arguments) == 1
-    V = arguments[0].ufl_function_space()
-
+    # Compile expression
     num_points = q_points.shape[0]
     compiled_expr = dolfinx.fem.Expression(expr, q_points)
-    mesh = Q.mesh
-    tdim = mesh.topology.dim
-    num_cells = mesh.topology.index_map(tdim).size_local
-    #
+
     # (num_cells, num_points, num_dofs*bs, expr_value_size)
-    array_evaluated = compiled_expr.eval(mesh, np.arange(num_cells, dtype=np.int32))
+    array_evaluated = compiled_expr.eval(mesh, interpolation_entities)
     assert np.prod(Q.value_shape) == np.prod(expr.ufl_shape)
 
     # Get data as (num_cells*num_points,1, expr_shape, num_test_basis_functions*test_block_size)