msh_0.geometry.dofmap.links(i) returns vertex ids of the ith cell, e.g. it returns 3 integers for mesh with triangle elements. To use that with the following
cell_geom_nodes = msh_0.geometry.dofmap.links(i)
cell_geom_coords = msh_0.geometry.x[cell_geom_nodes]
In v0.8.0, you can do
cell_geom_nodes = msh_0.geometry.dofmap[i]
cell_geom_coords = msh_0.geometry.x[cell_geom_nodes]
I also want to add that there’s a unit test for interpolation that includes nonmatching meshes, e.g. at dolfinx/python/test/unit/fem/test_interpolation.py at main · FEniCS/dolfinx · GitHub
I am sorry, but at the moment I am still working with version 0.6.0 as well. I have built an entire set of tools for vibroacoustics with it, and my daily work relies on it. I will port all the scripts to the most recent version, but I cannot promise when I will do that.
You can also have a look at the following post, where strategy and functionality are used to create a interpolation matrix in parallel.
I managed to get this algorithm working in dolfinx v0.9.0
import numpy as np
import basix
import dolfinx
from dolfinx import geometry
from scipy.sparse import coo_matrix
def interpolation_matrix_nonmatching_meshes(V_0, V_1, sparse=False, tol=1e-8):
"""
Compute the interpolation matrix between nonmatching meshes.
Parameters
----------
V_0 : dolfinx.FunctionSpace
Source function space.
V_1 : dolfinx.FunctionSpace
Target function space.
sparse : bool, optional
If True, return a SciPy CSR sparse matrix. Otherwise return a dense NumPy array.
tol : float, optional
Tolerance for checking point containment in reference element.
Returns
-------
I : numpy.ndarray or scipy.sparse.csr_matrix
Interpolation matrix R such that u_1 = R u_0. Type depends on `sparse` flag.
Raises
------
ValueError
If any target DOF point does not lie in a source cell.
"""
msh_0 = V_0.mesh
x_0 = V_0.tabulate_dof_coordinates()
x_1 = V_1.tabulate_dof_coordinates()
# Broad-phase collision detection
bb_tree = geometry.bb_tree(msh_0, msh_0.topology.dim)
cell_candidates = geometry.compute_collisions_points(bb_tree, x_1)
colliding_cells = geometry.compute_colliding_cells(msh_0, cell_candidates, x_1)
cells, points, idxs, ref_pts = [], [], [], []
missing = []
# Choose the best colliding cell based on containment or minimal violation
for i, point in enumerate(x_1):
links = colliding_cells.links(i)
# Correct check for no links
if len(links) == 0:
missing.append(i)
continue
best_cell = None
best_ref = None
best_violation = np.inf
for cidx in links:
gdofs = msh_0.geometry.dofmap[cidx]
cell_geo = msh_0.geometry.x[gdofs]
ref_pt = msh_0.geometry.cmap.pull_back(point.reshape(1, -1), cell_geo)[0]
# Check containment (simplex): barycentric coords >= -tol and sum <= 1+tol
if np.all(ref_pt >= -tol) and np.sum(ref_pt) <= 1.0 + tol:
best_cell = cidx
best_ref = ref_pt
break
else:
# Violation metric: sum of negative or exceeding portions
violation = (np.sum(np.maximum(-ref_pt, 0)) +
np.sum(np.maximum(ref_pt - 1.0, 0)) +
max(np.sum(ref_pt) - 1.0, 0))
if violation < best_violation:
best_violation = violation
best_cell = cidx
best_ref = ref_pt
cells.append(best_cell)
points.append(point)
idxs.append(i)
ref_pts.append(best_ref)
if missing:
raise ValueError(f"Target DOF indices with no source cell: {missing}")
idxs = np.array(idxs, dtype=int)
points = np.array(points, dtype=np.float64)
cells = np.array(cells, dtype=int)
x_ref = np.array(ref_pts, dtype=np.float64)
# Create finite element for evaluation
ct = dolfinx.cpp.mesh.to_string(msh_0.topology.cell_type)
element = basix.create_element(
basix.finite_element.string_to_family("Lagrange", ct),
msh_0.basix_cell(),
V_0.ufl_element().degree,
basix.LagrangeVariant.equispaced
)
# Evaluate basis at reference points
B = element.tabulate(0, x_ref)[0, :, :, 0]
ndofs_loc = B.shape[1]
# Build sparse triplet lists
rows, cols, data = [], [], []
for k, gi in enumerate(idxs):
cell_dofs = V_0.dofmap.cell_dofs(cells[k])
for j in range(ndofs_loc):
val = B[k, j]
if val != 0.0:
rows.append(gi)
cols.append(cell_dofs[j])
data.append(val)
# Return in requested format
if sparse:
I_sparse = coo_matrix((data, (rows, cols)), shape=(len(x_1), len(x_0)))
return I_sparse.tocsr()
else:
I = np.zeros((len(x_1), len(x_0)), dtype=np.float64)
for r, c, v in zip(rows, cols, data):
I[r, c] = v
return I