Return value meaning of nmm_interpolation_data

dolfinx
1

What is the meaning of the 4 return values of create_nonmatching_meshes_interpolation_data?

+----------------------+---------------------------+---------------------------------+
| C++ pybind11 m.def() | determine_point_ownership | interpolate_nonmatching_meshes  |
+----------------------+---------------------------+---------------------------------+
| src_owner            | point_owners              | dest_ranks                      |
| dest_owner           | owned_recv_ranks          | src_ranks                       |
| dest_points          | owned_recv_points         | recv_points                     |
| dest_cells           | owned_recv_cells          | evaluation_cells                |
+----------------------+---------------------------+---------------------------------+

I found them in different places as tabulated above and I guess the meaning based on the context:
(1) src_owner: the remote process rank with those points whose values are evaluated on the current process.
(2) dest_owner: a process rank, no clue…
(3) dest_points: coordinates of those points received from remote processes in (1).
(4) dest_cells: the cell indices of what dest_points belong to on the current process.

The relation between dest_points is dest_cell is len(dest_points) = 3 * len(dest_cells). And I obseved that len(dest_owner) = len(dest_cells).

2

We recently improved the documentation of that function, please see if Wrap `create_nonmatching_meshes_interpolation_data` in Python by mscroggs · Pull Request #3039 · FEniCS/dolfinx · GitHub helps.

1 Like
3

I am still lost with the meaning of src_owner. It says src_owner: Ranks owning each point sent into ownership determination for current process.

Actually, my initial understanding for src_owner should be dest_owner.

4

Imagine that you have two processes, rank 0 and rank 1

Rank 0 has 5 points, p0, p1, p2, p3, p4

Rank 1 has 2 points, q0, q1

You want to figure out which of the two processes owns the collection of points (p, q).
Let’s say Rank 0, owns p2, p4, q1, while Rank 1 owns p0, p1, p3, q0

Then src_owner on rank 0 is [1,1,0,1,0], while src_owner on rank 1 [1, 0].

When non-matching interpolation data is called, the aforementioned points are sent to a process (dest_points) and they know where they are coming from (dest_source).

5

A follow-up question: How can I associate the point with DOF of a FunctionSpace?

To brief the question: I do interpolation for two non-matching meshes in a serial run. The point info that create_nonmatching_meshes_interpolation_data returns is actually the coordinates of DOFs. I observed that the DOF points are retrieved by cells, which means a DOF on a vertex might be counted serveral times if it is shared by surrounding cells.

The following figure and code show what I do: Interpolate FunctionSpace Q1 on mesh1 into FunctionSpace Q2 on mesh2. I wonder how to get the DOF index of Q2 from the given point info?

import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh, cpp

# Mesh
n1 = 2
n2 = 1
mesh1 = mesh.create_unit_square(MPI.COMM_WORLD, n1, n1, cell_type=cpp.mesh.CellType.quadrilateral)
mesh2 = mesh.create_unit_square(MPI.COMM_WORLD, n2, n2, cell_type=cpp.mesh.CellType.triangle)
mesh2.geometry.x[:, : ] *= 0.5
mesh2.geometry.x[:, :2] += 0.25

# Function spaces
Q1 = fem.FunctionSpace(mesh1, ("Lagrange", 1))
Q2 = fem.FunctionSpace(mesh2, ("Lagrange", 1))

# Interpolation data
nms = fem.create_nonmatching_meshes_interpolation_data(
        Q2.mesh._cpp_object,
        Q2.element,
        Q1.mesh._cpp_object, 
        padding=1e-6)

point_owners, owned_recv_ranks, owned_recv_points, owned_recv_cells = nms

# Reshaped owned_recv_points -- rs_owned_recv_points
rs_owned_recv_points = np.reshape(owned_recv_points, (-1, 3))

# Geometry
geometry1 = mesh1.geometry
geometry2 = mesh2.geometry

# Geometry connectivity: cell
cells_1 = geometry1.dofmap
cells_2 = geometry2.dofmap

# Exam the point-cell relations
for i, cell in enumerate(owned_recv_cells):
    point = rs_owned_recv_points[i]
    cell_coords = geometry1.x[geometry1.dofmap[cell]]
    print(f"{i}th point of owned_recv_points with coordinate {point}" 
          f" falls in the {cell}th cell of mesh1 with coords\n{cell_coords}\n")

point2dof
point2dof1795×740 48.4 KB

6

Not all function spaces has a single point associated with a degree of freedom (they might have many, as a dof might be defined as integral moments).

It is unclear to me why this is an issue for you, as interpolation from one mesh to another is done on a cell to cell basis, i.e. for a given cell, get the interpolation points on the reference elements, push them forward to physical space, find an owner for them, evaluate the function and then return the values.

7

Try to concise my question: how to assign the evaluated point values back to corresponding DOFs efficiently?

# Build point to dof mapping, O(n^2)
point_to_dof = np.zeros(len(owned_recv_cells), dtype=np.int32)
for i, point in enumerate(rs_owned_recv_points):
    for j, dof_coordinates in enumerate(Q2.tabulate_dof_coordinates()):
        if np.linalg.norm(point - dof_coordinates) < 1e-6:
            point_to_dof[i] = j

The following code demonstrates my attempt: I interpolate u1 into u2 with create_nonmatching_meshes_interpolation_data. u2 is built on a mesh with 2 triangles and it has 4 DOFs. The number of point evaulated is 6. I build a point_to_dof mapping by comparing point and the DOF coordinates as shown above. However, this causes efficiency downgrade when the scale is large. Is there a better way of doing so or any advice to improve?

MWE:

import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh, cpp

# Mesh
mesh1 = mesh.create_unit_square(MPI.COMM_WORLD, 2, 2, cell_type=cpp.mesh.CellType.quadrilateral)
mesh2 = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, cell_type=cpp.mesh.CellType.triangle)
mesh2.geometry.x[:, : ] *= 0.5
mesh2.geometry.x[:, :2] += 0.25

# Function spaces and functions
Q1 = fem.FunctionSpace(mesh1, ("Lagrange", 1))
Q2 = fem.FunctionSpace(mesh2, ("Lagrange", 1))
u1 = fem.Function(Q1)
u2 = fem.Function(Q2)
u2_true = fem.Function(Q2, name="validation")
u1.interpolate(lambda x: np.sqrt(x[0]*x[0] + x[1]*x[1]))

# Interpolation data
nms = fem.create_nonmatching_meshes_interpolation_data(
        Q2.mesh._cpp_object,
        Q2.element,
        Q1.mesh._cpp_object, 
        padding=1e-6)

point_owners, owned_recv_ranks, owned_recv_points, owned_recv_cells = nms

# Reshaped owned_recv_points -- rs_owned_recv_points
rs_owned_recv_points = np.reshape(owned_recv_points, (-1, 3))

# Build point to dof mapping, O(n^2)
point_to_dof = np.zeros(len(owned_recv_cells), dtype=np.int32)
for i, point in enumerate(rs_owned_recv_points):
    for j, dof_coordinates in enumerate(Q2.tabulate_dof_coordinates()):
        if np.linalg.norm(point - dof_coordinates) < 1e-6:
            point_to_dof[i] = j

# `point` evaluation
evaluated_points = np.zeros(len(owned_recv_cells))
num = len(rs_owned_recv_points)
for i, point, cell in zip(range(num), rs_owned_recv_points, owned_recv_cells):
    evaluated_points[i] = u1.eval(point, cell)[0]

# Assign `point` value back to corresponding DOF
for i, evaluated_point in enumerate(evaluated_points):
    u2.x.array[point_to_dof[i]] = evaluated_point

# Validation
u2_true.interpolate(u1, nmm_interpolation_data=nms)
assert(np.allclose(u2.x.array, u2_true.x.array))
8

As I already stated, there isn’t always a 1-1 correspondence between the points and degrees of freedom. They have to be multiplied by an interpolation matrix for a large set of finite elements, as seen in:

the points are ordered in the following way: (number_of_cells, number_of_interpolation_points).

  1. Take each cell, access its dofmap to get the dofs of interest
  2. Get the interpolation points for the specific cell, (which in your case are those received from rs_owned_points
  3. Pull coefficients in current cell back, see section 8.4.3 of: DOLFINx: The next generation FEniCS problem solving environment
  4. Evaluate basis functions at given points, multiply pulled back dofs by interpolation matrix to get the dof values for given cell.
  5. Insert data in global matrix.

Please note that you could also build a vectorized distance matrix for your double for-loop, i.e. How to represent non-uniform first and second type boundary conditions? - #6 by dokken

9

I still want to check my understanding on the first return value point_owners. Please correct me if I am wrong.

It is a list of int which infers the process ranks to evaluate each DOF point that the current process owns.
Its length may be greater than the number of DOFs that the current process owns as the DOFs are gathered from each cell and piled up.
Then the length of point_owners should be equal to the following expression.

len(point_owners) = cell_map_total_num * function_space.dofmap.dof_layout.num_dofs

And the corresponding DOF IDs of the list point_owners should be retrived as follows.

point = [dof for i in range(cell_map_total_num) for dof in function_space.dofmap.cell_dofs(i)]

The full MWE:

# QnA_point_owners.py
# mpirun -n 2 python QnA_point_owners.py

from mpi4py import MPI
from dolfinx import fem, mesh, cpp

rank = MPI.COMM_WORLD.rank

# Mesh
n1, n2 = 5, 5
meshf = mesh.create_unit_square(MPI.COMM_WORLD, n1, n1, cell_type=cpp.mesh.CellType.quadrilateral)
meshs = mesh.create_unit_square(MPI.COMM_WORLD, n2, n2, cell_type=cpp.mesh.CellType.triangle)
meshs.geometry.x[:, : ] *= 0.5
meshs.geometry.x[:, :2] += 0.25

# Function spaces and functions
order = 1
Qf = fem.FunctionSpace(meshf, ("Lagrange", order))
Qs = fem.FunctionSpace(meshs, ("Lagrange", order))

# Interpolation data
point_owners, owned_recv_ranks, owned_recv_points, owned_recv_cells = (
    fem.create_nonmatching_meshes_interpolation_data(
        Qs.mesh._cpp_object,
        Qs.element,
        Qf.mesh._cpp_object, 
        padding=1e-6)
    )

topology_s = Qs.mesh.topology
tdim_s = Qs.mesh.topology.dim
cell_map_s = topology_s.index_map(tdim_s)
cell_total_num_s = cell_map_s.size_local + cell_map_s.num_ghosts

assert len(point_owners) == cell_total_num_s * Qs.dofmap.dof_layout.num_dofs
points = [dof for i in range(cell_total_num_s) for dof in Qs.dofmap.cell_dofs(i)]  

print(f"rank\n{rank}\n"
      f"len(point_owners)\n{len(point_owners)}\n"
      f"points:\n{points}\n")
10

Partially true. For a function space where each functional is a point evaluation (like Lagrange), this is true. However, for more exotic spaces, like Nedelec or RT, where dofs are defined through integral moments, a single dof might require multiple points (to create an approximation of this integral).

1 Like