dokken
March 7, 2024, 11:55am
8
chenyongxin:
# 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
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:
for (std::size_t c = 0; c < cells.size(); ++c)
{
const std::int32_t cell = cells[c];
std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
for (int k = 0; k < element_bs; ++k)
{
// num_scalar_dofs is the number of interpolation points per
// cell in this case (interpolation matrix is identity)
std::copy_n(std::next(f.begin(), k * f_shape1 + c * num_scalar_dofs),
num_scalar_dofs, _coeffs.begin());
apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
for (int i = 0; i < num_scalar_dofs; ++i)
{
const int dof = i * element_bs + k;
std::div_t pos = std::div(dof, dofmap_bs);
coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
}
}
}
}
the points are ordered in the following way: (number_of_cells, number_of_interpolation_points).
Take each cell, access its dofmap to get the dofs of interest
Get the interpolation points for the specific cell, (which in your case are those received from rs_owned_points
Pull coefficients in current cell back, see section 8.4.3 of: DOLFINx: The next generation FEniCS problem solving environment
Evaluate basis functions at given points, multiply pulled back dofs by interpolation matrix to get the dof values for given cell.
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