DOLFINx
DOLFINx C++ interface
Loading...
Searching...
No Matches
interpolate.h
1// Copyright (C) 2020-2022 Garth N. Wells, Igor A. Baratta, Massimiliano Leoni
2// and Jørgen S.Dokken
3//
4// This file is part of DOLFINx (https://www.fenicsproject.org)
5//
6// SPDX-License-Identifier: LGPL-3.0-or-later
7
8#pragma once
9
10#include "CoordinateElement.h"
11#include "DofMap.h"
12#include "FiniteElement.h"
13#include "FunctionSpace.h"
14#include <basix/mdspan.hpp>
15#include <concepts>
16#include <dolfinx/common/IndexMap.h>
17#include <dolfinx/common/types.h>
18#include <dolfinx/geometry/BoundingBoxTree.h>
19#include <dolfinx/geometry/utils.h>
20#include <dolfinx/mesh/Mesh.h>
21#include <functional>
22#include <numeric>
23#include <span>
24#include <vector>
25
26namespace dolfinx::fem
27{
28template <dolfinx::scalar T, std::floating_point U>
29class Function;
30
31template <typename T>
32concept MDSpan = requires(T x, std::size_t idx) {
33 x(idx, idx);
34 {
35 x.extent(0)
36 } -> std::integral;
37
38 {
39 x.extent(1)
40 } -> std::integral;
41};
42
53template <std::floating_point T>
54std::vector<T> interpolation_coords(const fem::FiniteElement<T>& element,
55 const mesh::Geometry<T>& geometry,
56 std::span<const std::int32_t> cells)
57{
58 // Get geometry data and the element coordinate map
59 const std::size_t gdim = geometry.dim();
60 auto x_dofmap = geometry.dofmap();
61 std::span<const T> x_g = geometry.x();
62
63 const CoordinateElement<T>& cmap = geometry.cmap();
64 const std::size_t num_dofs_g = cmap.dim();
65
66 // Get the interpolation points on the reference cells
67 const auto [X, Xshape] = element.interpolation_points();
68
69 // Evaluate coordinate element basis at reference points
70 std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(0, Xshape[0]);
71 std::vector<T> phi_b(
72 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
73 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
74 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>
75 phi_full(phi_b.data(), phi_shape);
76 cmap.tabulate(0, X, Xshape, phi_b);
77 auto phi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
78 phi_full, 0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
79 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
80
81 // Push reference coordinates (X) forward to the physical coordinates
82 // (x) for each cell
83 std::vector<T> coordinate_dofs(num_dofs_g * gdim, 0);
84 std::vector<T> x(3 * (cells.size() * Xshape[0]), 0);
85 for (std::size_t c = 0; c < cells.size(); ++c)
86 {
87 // Get geometry data for current cell
88 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
89 x_dofmap, cells[c], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
90 for (std::size_t i = 0; i < x_dofs.size(); ++i)
91 {
92 std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), gdim,
93 std::next(coordinate_dofs.begin(), i * gdim));
94 }
95
96 // Push forward coordinates (X -> x)
97 for (std::size_t p = 0; p < Xshape[0]; ++p)
98 {
99 for (std::size_t j = 0; j < gdim; ++j)
100 {
101 T acc = 0;
102 for (std::size_t k = 0; k < num_dofs_g; ++k)
103 acc += phi(p, k) * coordinate_dofs[k * gdim + j];
104 x[j * (cells.size() * Xshape[0]) + c * Xshape[0] + p] = acc;
105 }
106 }
107 }
108
109 return x;
110}
111
126template <dolfinx::scalar T, std::floating_point U>
127void interpolate(Function<T, U>& u, std::span<const T> f,
128 std::array<std::size_t, 2> fshape,
129 std::span<const std::int32_t> cells);
130
131namespace impl
132{
134template <typename T, std::size_t D>
135using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
136 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, D>>;
137
157template <dolfinx::scalar T>
158void scatter_values(
159 MPI_Comm comm, std::span<const std::int32_t> src_ranks,
160 std::span<const std::int32_t> dest_ranks,
161 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
162 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
163 send_values,
164 std::span<T> recv_values)
165{
166 const std::size_t block_size = send_values.extent(1);
167 assert(src_ranks.size() * block_size == send_values.size());
168 assert(recv_values.size() == dest_ranks.size() * block_size);
169
170 // Build unique set of the sorted src_ranks
171 std::vector<std::int32_t> out_ranks(src_ranks.size());
172 out_ranks.assign(src_ranks.begin(), src_ranks.end());
173 out_ranks.erase(std::unique(out_ranks.begin(), out_ranks.end()),
174 out_ranks.end());
175 out_ranks.reserve(out_ranks.size() + 1);
176
177 // Remove negative entries from dest_ranks
178 std::vector<std::int32_t> in_ranks;
179 in_ranks.reserve(dest_ranks.size());
180 std::copy_if(dest_ranks.begin(), dest_ranks.end(),
181 std::back_inserter(in_ranks),
182 [](auto rank) { return rank >= 0; });
183
184 // Create unique set of sorted in-ranks
185 std::sort(in_ranks.begin(), in_ranks.end());
186 in_ranks.erase(std::unique(in_ranks.begin(), in_ranks.end()), in_ranks.end());
187 in_ranks.reserve(in_ranks.size() + 1);
188
189 // Create neighborhood communicator
190 MPI_Comm reverse_comm;
191 MPI_Dist_graph_create_adjacent(
192 comm, in_ranks.size(), in_ranks.data(), MPI_UNWEIGHTED, out_ranks.size(),
193 out_ranks.data(), MPI_UNWEIGHTED, MPI_INFO_NULL, false, &reverse_comm);
194
195 std::vector<std::int32_t> comm_to_output;
196 std::vector<std::int32_t> recv_sizes(in_ranks.size());
197 recv_sizes.reserve(1);
198 std::vector<std::int32_t> recv_offsets(in_ranks.size() + 1, 0);
199 {
200 // Build map from parent to neighborhood communicator ranks
201 std::map<std::int32_t, std::int32_t> rank_to_neighbor;
202 for (std::size_t i = 0; i < in_ranks.size(); i++)
203 rank_to_neighbor[in_ranks[i]] = i;
204
205 // Compute receive sizes
206 std::for_each(
207 dest_ranks.begin(), dest_ranks.end(),
208 [&dest_ranks, &rank_to_neighbor, &recv_sizes, block_size](auto rank)
209 {
210 if (rank >= 0)
211 {
212 const int neighbor = rank_to_neighbor[rank];
213 recv_sizes[neighbor] += block_size;
214 }
215 });
216
217 // Compute receiving offsets
218 std::partial_sum(recv_sizes.begin(), recv_sizes.end(),
219 std::next(recv_offsets.begin(), 1));
220
221 // Compute map from receiving values to position in recv_values
222 comm_to_output.resize(recv_offsets.back() / block_size);
223 std::vector<std::int32_t> recv_counter(recv_sizes.size(), 0);
224 for (std::size_t i = 0; i < dest_ranks.size(); ++i)
225 {
226 if (const std::int32_t rank = dest_ranks[i]; rank >= 0)
227 {
228 const int neighbor = rank_to_neighbor[rank];
229 int insert_pos = recv_offsets[neighbor] + recv_counter[neighbor];
230 comm_to_output[insert_pos / block_size] = i * block_size;
231 recv_counter[neighbor] += block_size;
232 }
233 }
234 }
235
236 std::vector<std::int32_t> send_sizes(out_ranks.size());
237 send_sizes.reserve(1);
238 {
239 // Compute map from parent mpi rank to neigbor rank for outgoing data
240 std::map<std::int32_t, std::int32_t> rank_to_neighbor;
241 for (std::size_t i = 0; i < out_ranks.size(); i++)
242 rank_to_neighbor[out_ranks[i]] = i;
243
244 // Compute send sizes
245 std::for_each(src_ranks.begin(), src_ranks.end(),
246 [&rank_to_neighbor, &send_sizes, block_size](auto rank)
247 { send_sizes[rank_to_neighbor[rank]] += block_size; });
248 }
249
250 // Compute sending offsets
251 std::vector<std::int32_t> send_offsets(send_sizes.size() + 1, 0);
252 std::partial_sum(send_sizes.begin(), send_sizes.end(),
253 std::next(send_offsets.begin(), 1));
254
255 // Send values to dest ranks
256 std::vector<T> values(recv_offsets.back());
257 values.reserve(1);
258 MPI_Neighbor_alltoallv(send_values.data_handle(), send_sizes.data(),
259 send_offsets.data(), dolfinx::MPI::mpi_type<T>(),
260 values.data(), recv_sizes.data(), recv_offsets.data(),
261 dolfinx::MPI::mpi_type<T>(), reverse_comm);
262 MPI_Comm_free(&reverse_comm);
263
264 // Insert values received from neighborhood communicator in output span
265 std::fill(recv_values.begin(), recv_values.end(), T(0));
266 for (std::size_t i = 0; i < comm_to_output.size(); i++)
267 {
268 auto vals = std::next(recv_values.begin(), comm_to_output[i]);
269 auto vals_from = std::next(values.begin(), i * block_size);
270 std::copy_n(vals_from, block_size, vals);
271 }
272};
273
282template <MDSpan U, MDSpan V, dolfinx::scalar T>
283void interpolation_apply(U&& Pi, V&& data, std::span<T> coeffs, int bs)
284{
285 using X = typename dolfinx::scalar_value_type_t<T>;
286
287 // Compute coefficients = Pi * x (matrix-vector multiply)
288 if (bs == 1)
289 {
290 assert(data.extent(0) * data.extent(1) == Pi.extent(1));
291 for (std::size_t i = 0; i < Pi.extent(0); ++i)
292 {
293 coeffs[i] = 0.0;
294 for (std::size_t k = 0; k < data.extent(1); ++k)
295 for (std::size_t j = 0; j < data.extent(0); ++j)
296 coeffs[i]
297 += static_cast<X>(Pi(i, k * data.extent(0) + j)) * data(j, k);
298 }
299 }
300 else
301 {
302 const std::size_t cols = Pi.extent(1);
303 assert(data.extent(0) == Pi.extent(1));
304 assert(data.extent(1) == bs);
305 for (int k = 0; k < bs; ++k)
306 {
307 for (std::size_t i = 0; i < Pi.extent(0); ++i)
308 {
309 T acc = 0;
310 for (std::size_t j = 0; j < cols; ++j)
311 acc += static_cast<X>(Pi(i, j)) * data(j, k);
312 coeffs[bs * i + k] = acc;
313 }
314 }
315 }
316}
317
328template <dolfinx::scalar T, std::floating_point U>
329void interpolate_same_map(Function<T, U>& u1, const Function<T, U>& u0,
330 std::span<const std::int32_t> cells)
331{
332 auto V0 = u0.function_space();
333 assert(V0);
334 auto V1 = u1.function_space();
335 assert(V1);
336 auto mesh = V0->mesh();
337 assert(mesh);
338
339 auto element0 = V0->element();
340 assert(element0);
341 auto element1 = V1->element();
342 assert(element1);
343
344 const int tdim = mesh->topology()->dim();
345 auto map = mesh->topology()->index_map(tdim);
346 assert(map);
347 std::span<T> u1_array = u1.x()->mutable_array();
348 std::span<const T> u0_array = u0.x()->array();
349
350 std::span<const std::uint32_t> cell_info;
351 if (element1->needs_dof_transformations()
352 or element0->needs_dof_transformations())
353 {
354 mesh->topology_mutable()->create_entity_permutations();
355 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
356 }
357
358 // Get dofmaps
359 auto dofmap1 = V1->dofmap();
360 auto dofmap0 = V0->dofmap();
361
362 // Get block sizes and dof transformation operators
363 const int bs1 = dofmap1->bs();
364 const int bs0 = dofmap0->bs();
365 auto apply_dof_transformation
366 = element0->template get_pre_dof_transformation_function<T>(
367 FiniteElement<U>::doftransform::transpose, false);
368 auto apply_inverse_dof_transform
369 = element1->template get_pre_dof_transformation_function<T>(
370 FiniteElement<U>::doftransform::inverse_transpose, false);
371
372 // Create working array
373 std::vector<T> local0(element0->space_dimension());
374 std::vector<T> local1(element1->space_dimension());
375
376 // Create interpolation operator
377 const auto [i_m, im_shape]
378 = element1->create_interpolation_operator(*element0);
379
380 // Iterate over mesh and interpolate on each cell
381 using X = typename dolfinx::scalar_value_type_t<T>;
382 for (auto c : cells)
383 {
384 std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
385 for (std::size_t i = 0; i < dofs0.size(); ++i)
386 for (int k = 0; k < bs0; ++k)
387 local0[bs0 * i + k] = u0_array[bs0 * dofs0[i] + k];
388
389 apply_dof_transformation(local0, cell_info, c, 1);
390
391 // FIXME: Get compile-time ranges from Basix
392 // Apply interpolation operator
393 std::fill(local1.begin(), local1.end(), 0);
394 for (std::size_t i = 0; i < im_shape[0]; ++i)
395 for (std::size_t j = 0; j < im_shape[1]; ++j)
396 local1[i] += static_cast<X>(i_m[im_shape[1] * i + j]) * local0[j];
397
398 apply_inverse_dof_transform(local1, cell_info, c, 1);
399 std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
400 for (std::size_t i = 0; i < dofs1.size(); ++i)
401 for (int k = 0; k < bs1; ++k)
402 u1_array[bs1 * dofs1[i] + k] = local1[bs1 * i + k];
403 }
404}
405
416template <dolfinx::scalar T, std::floating_point U>
417void interpolate_nonmatching_maps(Function<T, U>& u1, const Function<T, U>& u0,
418 std::span<const std::int32_t> cells)
419{
420 // Get mesh
421 auto V0 = u0.function_space();
422 assert(V0);
423 auto mesh = V0->mesh();
424 assert(mesh);
425
426 // Mesh dims
427 const int tdim = mesh->topology()->dim();
428 const int gdim = mesh->geometry().dim();
429
430 // Get elements
431 auto V1 = u1.function_space();
432 assert(V1);
433 auto element0 = V0->element();
434 assert(element0);
435 auto element1 = V1->element();
436 assert(element1);
437
438 std::span<const std::uint32_t> cell_info;
439 if (element1->needs_dof_transformations()
440 or element0->needs_dof_transformations())
441 {
442 mesh->topology_mutable()->create_entity_permutations();
443 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
444 }
445
446 // Get dofmaps
447 auto dofmap0 = V0->dofmap();
448 auto dofmap1 = V1->dofmap();
449
450 const auto [X, Xshape] = element1->interpolation_points();
451
452 // Get block sizes and dof transformation operators
453 const int bs0 = element0->block_size();
454 const int bs1 = element1->block_size();
455 auto apply_dof_transformation0
456 = element0->template get_pre_dof_transformation_function<U>(
457 FiniteElement<U>::doftransform::standard, false);
458 auto apply_inverse_dof_transform1
459 = element1->template get_pre_dof_transformation_function<T>(
460 FiniteElement<U>::doftransform::inverse_transpose, false);
461
462 // Get sizes of elements
463 const std::size_t dim0 = element0->space_dimension() / bs0;
464 const std::size_t value_size_ref0 = element0->reference_value_size() / bs0;
465 const std::size_t value_size0 = V0->value_size() / bs0;
466
467 const CoordinateElement<U>& cmap = mesh->geometry().cmap();
468 auto x_dofmap = mesh->geometry().dofmap();
469 const std::size_t num_dofs_g = cmap.dim();
470 std::span<const U> x_g = mesh->geometry().x();
471
472 // Evaluate coordinate map basis at reference interpolation points
473 const std::array<std::size_t, 4> phi_shape
474 = cmap.tabulate_shape(1, Xshape[0]);
475 std::vector<U> phi_b(
476 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
477 mdspan_t<const U, 4> phi(phi_b.data(), phi_shape);
478 cmap.tabulate(1, X, Xshape, phi_b);
479
480 // Evaluate v basis functions at reference interpolation points
481 const auto [_basis_derivatives_reference0, b0shape]
482 = element0->tabulate(X, Xshape, 0);
483 mdspan_t<const U, 4> basis_derivatives_reference0(
484 _basis_derivatives_reference0.data(), b0shape);
485
486 // Create working arrays
487 std::vector<T> local1(element1->space_dimension());
488 std::vector<T> coeffs0(element0->space_dimension());
489
490 std::vector<U> basis0_b(Xshape[0] * dim0 * value_size0);
491 impl::mdspan_t<U, 3> basis0(basis0_b.data(), Xshape[0], dim0, value_size0);
492
493 std::vector<U> basis_reference0_b(Xshape[0] * dim0 * value_size_ref0);
494 impl::mdspan_t<U, 3> basis_reference0(basis_reference0_b.data(), Xshape[0],
495 dim0, value_size_ref0);
496
497 std::vector<T> values0_b(Xshape[0] * 1 * V1->value_size());
498 impl::mdspan_t<T, 3> values0(values0_b.data(), Xshape[0], 1,
499 V1->value_size());
500
501 std::vector<T> mapped_values_b(Xshape[0] * 1 * V1->value_size());
502 impl::mdspan_t<T, 3> mapped_values0(mapped_values_b.data(), Xshape[0], 1,
503 V1->value_size());
504
505 std::vector<U> coord_dofs_b(num_dofs_g * gdim);
506 impl::mdspan_t<U, 2> coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
507
508 std::vector<U> J_b(Xshape[0] * gdim * tdim);
509 impl::mdspan_t<U, 3> J(J_b.data(), Xshape[0], gdim, tdim);
510 std::vector<U> K_b(Xshape[0] * tdim * gdim);
511 impl::mdspan_t<U, 3> K(K_b.data(), Xshape[0], tdim, gdim);
512 std::vector<U> detJ(Xshape[0]);
513 std::vector<U> det_scratch(2 * gdim * tdim);
514
515 // Get interpolation operator
516 const auto [_Pi_1, pi_shape] = element1->interpolation_operator();
517 impl::mdspan_t<const U, 2> Pi_1(_Pi_1.data(), pi_shape);
518
519 using u_t = impl::mdspan_t<U, 2>;
520 using U_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
521 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
522 using J_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
523 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
524 using K_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
525 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
526 auto push_forward_fn0
527 = element0->basix_element().template map_fn<u_t, U_t, J_t, K_t>();
528
529 using v_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
530 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
531 using V_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
532 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
533 auto pull_back_fn1
534 = element1->basix_element().template map_fn<V_t, v_t, K_t, J_t>();
535
536 // Iterate over mesh and interpolate on each cell
537 std::span<const T> array0 = u0.x()->array();
538 std::span<T> array1 = u1.x()->mutable_array();
539 for (auto c : cells)
540 {
541 // Get cell geometry (coordinate dofs)
542 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
543 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
544 for (std::size_t i = 0; i < num_dofs_g; ++i)
545 {
546 const int pos = 3 * x_dofs[i];
547 for (std::size_t j = 0; j < gdim; ++j)
548 coord_dofs(i, j) = x_g[pos + j];
549 }
550
551 // Compute Jacobians and reference points for current cell
552 std::fill(J_b.begin(), J_b.end(), 0);
553 for (std::size_t p = 0; p < Xshape[0]; ++p)
554 {
555 auto dphi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
556 phi, std::pair(1, tdim + 1), p,
557 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
558
559 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
560 J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
561 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
562 cmap.compute_jacobian(dphi, coord_dofs, _J);
563 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
564 K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
565 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
566 cmap.compute_jacobian_inverse(_J, _K);
567 detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
568 }
569
570 // Copy evaluated basis on reference, apply DOF transformations, and
571 // push forward to physical element
572 for (std::size_t k0 = 0; k0 < basis_reference0.extent(0); ++k0)
573 for (std::size_t k1 = 0; k1 < basis_reference0.extent(1); ++k1)
574 for (std::size_t k2 = 0; k2 < basis_reference0.extent(2); ++k2)
575 basis_reference0(k0, k1, k2)
576 = basis_derivatives_reference0(0, k0, k1, k2);
577
578 for (std::size_t p = 0; p < Xshape[0]; ++p)
579 {
580 apply_dof_transformation0(
581 std::span(basis_reference0_b.data() + p * dim0 * value_size_ref0,
582 dim0 * value_size_ref0),
583 cell_info, c, value_size_ref0);
584 }
585
586 for (std::size_t i = 0; i < basis0.extent(0); ++i)
587 {
588 auto _u = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
589 basis0, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
590 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
591 auto _U = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
592 basis_reference0, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
593 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
594 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
595 K, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
596 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
597 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
598 J, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
599 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
600 push_forward_fn0(_u, _U, _J, detJ[i], _K);
601 }
602
603 // Copy expansion coefficients for v into local array
604 const int dof_bs0 = dofmap0->bs();
605 std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
606 for (std::size_t i = 0; i < dofs0.size(); ++i)
607 for (int k = 0; k < dof_bs0; ++k)
608 coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k];
609
610 // Evaluate v at the interpolation points (physical space values)
611 using X = typename dolfinx::scalar_value_type_t<T>;
612 for (std::size_t p = 0; p < Xshape[0]; ++p)
613 {
614 for (int k = 0; k < bs0; ++k)
615 {
616 for (std::size_t j = 0; j < value_size0; ++j)
617 {
618 T acc = 0;
619 for (std::size_t i = 0; i < dim0; ++i)
620 acc += coeffs0[bs0 * i + k] * static_cast<X>(basis0(p, i, j));
621 values0(p, 0, j * bs0 + k) = acc;
622 }
623 }
624 }
625
626 // Pull back the physical values to the u reference
627 for (std::size_t i = 0; i < values0.extent(0); ++i)
628 {
629 auto _u = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
630 values0, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
631 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
632 auto _U = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
633 mapped_values0, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
634 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
635 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
636 K, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
637 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
638 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
639 J, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
640 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
641 pull_back_fn1(_U, _u, _K, 1.0 / detJ[i], _J);
642 }
643
644 auto values = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
645 mapped_values0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0,
646 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
647 interpolation_apply(Pi_1, values, std::span(local1), bs1);
648 apply_inverse_dof_transform1(local1, cell_info, c, 1);
649
650 // Copy local coefficients to the correct position in u dof array
651 const int dof_bs1 = dofmap1->bs();
652 std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
653 for (std::size_t i = 0; i < dofs1.size(); ++i)
654 for (int k = 0; k < dof_bs1; ++k)
655 array1[dof_bs1 * dofs1[i] + k] = local1[dof_bs1 * i + k];
656 }
657}
658
659template <dolfinx::scalar T, std::floating_point U>
660void interpolate_nonmatching_meshes(
661 Function<T, U>& u, const Function<T, U>& v,
662 std::span<const std::int32_t> cells,
663 const std::tuple<std::span<const std::int32_t>,
664 std::span<const std::int32_t>, std::span<const U>,
665 std::span<const std::int32_t>>& nmm_interpolation_data)
666{
667 auto mesh = u.function_space()->mesh();
668 assert(mesh);
669 MPI_Comm comm = mesh->comm();
670
671 {
672 auto mesh_v = v.function_space()->mesh();
673 assert(mesh_v);
674 int result;
675 MPI_Comm_compare(comm, mesh_v->comm(), &result);
676 if (result == MPI_UNEQUAL)
677 {
678 throw std::runtime_error("Interpolation on different meshes is only "
679 "supported with the same communicator.");
680 }
681 }
682
683 assert(mesh->topology());
684 auto cell_map = mesh->topology()->index_map(mesh->topology()->dim());
685 assert(cell_map);
686
687 auto element_u = u.function_space()->element();
688 assert(element_u);
689 const std::size_t value_size = u.function_space()->value_size();
690
691 const auto& [dest_ranks, src_ranks, recv_points, evaluation_cells]
692 = nmm_interpolation_data;
693
694 // Evaluate the interpolating function where possible
695 std::vector<T> send_values(recv_points.size() / 3 * value_size);
696 v.eval(recv_points, {recv_points.size() / 3, (std::size_t)3},
697 evaluation_cells, send_values, {recv_points.size() / 3, value_size});
698
699 // Send values back to owning process
700 std::vector<T> values_b(dest_ranks.size() * value_size);
701 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
702 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
703 _send_values(send_values.data(), src_ranks.size(), value_size);
704 impl::scatter_values(comm, src_ranks, dest_ranks, _send_values,
705 std::span(values_b));
706
707 // Transpose received data
708 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
709 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
710 values(values_b.data(), dest_ranks.size(), value_size);
711 std::vector<T> valuesT_b(value_size * dest_ranks.size());
712 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
713 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
714 valuesT(valuesT_b.data(), value_size, dest_ranks.size());
715 for (std::size_t i = 0; i < values.extent(0); ++i)
716 for (std::size_t j = 0; j < values.extent(1); ++j)
717 valuesT(j, i) = values(i, j);
718
719 // Call local interpolation operator
720 fem::interpolate<T>(u, valuesT_b, {valuesT.extent(0), valuesT.extent(1)},
721 cells);
722}
723//----------------------------------------------------------------------------
724} // namespace impl
725
726template <dolfinx::scalar T, std::floating_point U>
727void interpolate(Function<T, U>& u, std::span<const T> f,
728 std::array<std::size_t, 2> fshape,
729 std::span<const std::int32_t> cells)
730{
731 using cmdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
732 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
733 using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
734 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
735 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
736 U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
737 using mdspan3_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
738 U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>;
739 auto element = u.function_space()->element();
740 assert(element);
741 const int element_bs = element->block_size();
742 if (int num_sub = element->num_sub_elements();
743 num_sub > 0 and num_sub != element_bs)
744 {
745 throw std::runtime_error("Cannot directly interpolate a mixed space. "
746 "Interpolate into subspaces.");
747 }
748
749 // Get mesh
750 assert(u.function_space());
751 auto mesh = u.function_space()->mesh();
752 assert(mesh);
753
754 const int gdim = mesh->geometry().dim();
755 const int tdim = mesh->topology()->dim();
756
757 if (fshape[0] != (std::size_t)u.function_space()->value_size())
758 throw std::runtime_error("Interpolation data has the wrong shape/size.");
759
760 std::span<const std::uint32_t> cell_info;
761 if (element->needs_dof_transformations())
762 {
763 mesh->topology_mutable()->create_entity_permutations();
764 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
765 }
766
767 const std::size_t f_shape1 = f.size() / u.function_space()->value_size();
768 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
769 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
770 _f(f.data(), fshape);
771
772 // Get dofmap
773 const auto dofmap = u.function_space()->dofmap();
774 assert(dofmap);
775 const int dofmap_bs = dofmap->bs();
776
777 // Loop over cells and compute interpolation dofs
778 const int num_scalar_dofs = element->space_dimension() / element_bs;
779 const int value_size = u.function_space()->value_size() / element_bs;
780
781 std::span<T> coeffs = u.x()->mutable_array();
782 std::vector<T> _coeffs(num_scalar_dofs);
783
784 // This assumes that any element with an identity interpolation matrix
785 // is a point evaluation
786 if (element->map_ident() && element->interpolation_ident())
787 {
788 // Point evaluation element *and* the geometric map is the identity,
789 // e.g. not Piola mapped
790
791 auto apply_inv_transpose_dof_transformation
792 = element->template get_pre_dof_transformation_function<T>(
794
795 // Loop over cells
796 for (std::size_t c = 0; c < cells.size(); ++c)
797 {
798 const std::int32_t cell = cells[c];
799 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
800 for (int k = 0; k < element_bs; ++k)
801 {
802 // num_scalar_dofs is the number of interpolation points per
803 // cell in this case (interpolation matrix is identity)
804 std::copy_n(std::next(f.begin(), k * f_shape1 + c * num_scalar_dofs),
805 num_scalar_dofs, _coeffs.begin());
806 apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
807 for (int i = 0; i < num_scalar_dofs; ++i)
808 {
809 const int dof = i * element_bs + k;
810 std::div_t pos = std::div(dof, dofmap_bs);
811 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
812 }
813 }
814 }
815 }
816 else if (element->map_ident())
817 {
818 // Not a point evaluation, but the geometric map is the identity,
819 // e.g. not Piola mapped
820
821 const int element_vs = u.function_space()->value_size() / element_bs;
822
823 if (element_vs > 1 && element_bs > 1)
824 {
825 throw std::runtime_error(
826 "Interpolation into this element not supported.");
827 }
828
829 // Get interpolation operator
830 const auto [_Pi, pi_shape] = element->interpolation_operator();
831 cmdspan2_t Pi(_Pi.data(), pi_shape);
832 const std::size_t num_interp_points = Pi.extent(1);
833 assert(Pi.extent(0) == num_scalar_dofs);
834
835 auto apply_inv_transpose_dof_transformation
836 = element->template get_pre_dof_transformation_function<T>(
838
839 // Loop over cells
840 std::vector<T> ref_data_b(num_interp_points);
841 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
842 T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents<
843 std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 1>>
844 ref_data(ref_data_b.data(), num_interp_points, 1);
845 for (std::size_t c = 0; c < cells.size(); ++c)
846 {
847 const std::int32_t cell = cells[c];
848 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
849 for (int k = 0; k < element_bs; ++k)
850 {
851 for (int i = 0; i < element_vs; ++i)
852 {
853 std::copy_n(
854 std::next(f.begin(), (i + k) * f_shape1
855 + c * num_interp_points / element_vs),
856 num_interp_points / element_vs,
857 std::next(ref_data_b.begin(),
858 i * num_interp_points / element_vs));
859 }
860 impl::interpolation_apply(Pi, ref_data, std::span(_coeffs), 1);
861 apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
862 for (int i = 0; i < num_scalar_dofs; ++i)
863 {
864 const int dof = i * element_bs + k;
865 std::div_t pos = std::div(dof, dofmap_bs);
866 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
867 }
868 }
869 }
870 }
871 else
872 {
873 // Get the interpolation points on the reference cells
874 const auto [X, Xshape] = element->interpolation_points();
875 if (X.empty())
876 {
877 throw std::runtime_error(
878 "Interpolation into this space is not yet supported.");
879 }
880
881 if (_f.extent(1) != cells.size() * Xshape[0])
882 throw std::runtime_error("Interpolation data has the wrong shape.");
883
884 // Get coordinate map
885 const CoordinateElement<U>& cmap = mesh->geometry().cmap();
886
887 // Get geometry data
888 auto x_dofmap = mesh->geometry().dofmap();
889 const int num_dofs_g = cmap.dim();
890 std::span<const U> x_g = mesh->geometry().x();
891
892 // Create data structures for Jacobian info
893 std::vector<U> J_b(Xshape[0] * gdim * tdim);
894 mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
895 std::vector<U> K_b(Xshape[0] * tdim * gdim);
896 mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
897 std::vector<U> detJ(Xshape[0]);
898 std::vector<U> det_scratch(2 * gdim * tdim);
899
900 std::vector<U> coord_dofs_b(num_dofs_g * gdim);
901 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
902
903 std::vector<T> ref_data_b(Xshape[0] * 1 * value_size);
904 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
905 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
906 ref_data(ref_data_b.data(), Xshape[0], 1, value_size);
907
908 std::vector<T> _vals_b(Xshape[0] * 1 * value_size);
909 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
910 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
911 _vals(_vals_b.data(), Xshape[0], 1, value_size);
912
913 // Tabulate 1st derivative of shape functions at interpolation
914 // coords
915 std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(1, Xshape[0]);
916 std::vector<U> phi_b(
917 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
918 cmdspan4_t phi(phi_b.data(), phi_shape);
919 cmap.tabulate(1, X, Xshape, phi_b);
920 auto dphi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
921 phi, std::pair(1, tdim + 1),
922 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
923 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
924
925 const std::function<void(std::span<T>, std::span<const std::uint32_t>,
926 std::int32_t, int)>
927 apply_inverse_transpose_dof_transformation
928 = element->template get_pre_dof_transformation_function<T>(
930
931 // Get interpolation operator
932 const auto [_Pi, pi_shape] = element->interpolation_operator();
933 cmdspan2_t Pi(_Pi.data(), pi_shape);
934
935 using u_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
936 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
937 using U_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
938 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
939 using J_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
940 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
941 using K_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
942 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
943 auto pull_back_fn
944 = element->basix_element().template map_fn<U_t, u_t, J_t, K_t>();
945
946 for (std::size_t c = 0; c < cells.size(); ++c)
947 {
948 const std::int32_t cell = cells[c];
949 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
950 x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
951 for (int i = 0; i < num_dofs_g; ++i)
952 {
953 const int pos = 3 * x_dofs[i];
954 for (int j = 0; j < gdim; ++j)
955 coord_dofs(i, j) = x_g[pos + j];
956 }
957
958 // Compute J, detJ and K
959 std::fill(J_b.begin(), J_b.end(), 0);
960 for (std::size_t p = 0; p < Xshape[0]; ++p)
961 {
962 auto _dphi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
963 dphi, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, p,
964 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
965 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
966 J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
967 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
968 cmap.compute_jacobian(_dphi, coord_dofs, _J);
969 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
970 K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
971 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
972 cmap.compute_jacobian_inverse(_J, _K);
973 detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
974 }
975
976 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
977 for (int k = 0; k < element_bs; ++k)
978 {
979 // Extract computed expression values for element block k
980 for (int m = 0; m < value_size; ++m)
981 {
982 for (std::size_t k0 = 0; k0 < Xshape[0]; ++k0)
983 {
984 _vals(k0, 0, m)
985 = f[f_shape1 * (k * value_size + m) + c * Xshape[0] + k0];
986 }
987 }
988
989 // Get element degrees of freedom for block
990 for (std::size_t i = 0; i < Xshape[0]; ++i)
991 {
992 auto _u = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
993 _vals, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
994 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
995 auto _U = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
996 ref_data, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
997 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
998 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
999 K, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
1000 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
1001 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
1002 J, i, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
1003 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
1004 pull_back_fn(_U, _u, _K, 1.0 / detJ[i], _J);
1005 }
1006
1007 auto ref = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
1008 ref_data, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0,
1009 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
1010 impl::interpolation_apply(Pi, ref, std::span(_coeffs), element_bs);
1011 apply_inverse_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
1012
1013 // Copy interpolation dofs into coefficient vector
1014 assert(_coeffs.size() == num_scalar_dofs);
1015 for (int i = 0; i < num_scalar_dofs; ++i)
1016 {
1017 const int dof = i * element_bs + k;
1018 std::div_t pos = std::div(dof, dofmap_bs);
1019 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
1020 }
1021 }
1022 }
1023 }
1024}
1025
1043template <std::floating_point T>
1044std::tuple<std::vector<std::int32_t>, std::vector<std::int32_t>, std::vector<T>,
1045 std::vector<std::int32_t>>
1047 const mesh::Geometry<T>& geometry0, const FiniteElement<T>& element0,
1048 const mesh::Mesh<T>& mesh1, std::span<const std::int32_t> cells, T padding)
1049{
1050 // Collect all the points at which values are needed to define the
1051 // interpolating function
1052 const std::vector<T> coords
1053 = interpolation_coords(element0, geometry0, cells);
1054
1055 // Transpose interpolation coords
1056 std::vector<T> x(coords.size());
1057 std::size_t num_points = coords.size() / 3;
1058 for (std::size_t i = 0; i < num_points; ++i)
1059 for (std::size_t j = 0; j < 3; ++j)
1060 x[3 * i + j] = coords[i + j * num_points];
1061
1062 // Determine ownership of each point
1063 return geometry::determine_point_ownership<T>(mesh1, x, padding);
1064}
1065
1079template <std::floating_point T>
1080std::tuple<std::vector<std::int32_t>, std::vector<std::int32_t>, std::vector<T>,
1081 std::vector<std::int32_t>>
1083 const FiniteElement<T>& element0,
1084 const mesh::Mesh<T>& mesh1,
1085 T padding)
1086{
1087 int tdim = mesh0.topology()->dim();
1088 auto cell_map = mesh0.topology()->index_map(tdim);
1089 assert(cell_map);
1090 std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
1091 std::vector<std::int32_t> cells(num_cells, 0);
1092 std::iota(cells.begin(), cells.end(), 0);
1094 mesh0.geometry(), element0, mesh1, cells, padding);
1095}
1096
1104template <dolfinx::scalar T, std::floating_point U>
1106 Function<T, U>& u, const Function<T, U>& v,
1107 std::span<const std::int32_t> cells,
1108 const std::tuple<std::span<const std::int32_t>,
1109 std::span<const std::int32_t>, std::span<const U>,
1110 std::span<const std::int32_t>>& nmm_interpolation_data
1111 = {})
1112{
1113 assert(u.function_space());
1114 assert(v.function_space());
1115 auto mesh = u.function_space()->mesh();
1116 assert(mesh);
1117
1118 auto cell_map0 = mesh->topology()->index_map(mesh->topology()->dim());
1119 assert(cell_map0);
1120 std::size_t num_cells0 = cell_map0->size_local() + cell_map0->num_ghosts();
1121 if (u.function_space() == v.function_space() and cells.size() == num_cells0)
1122 {
1123 // Same function spaces and on whole mesh
1124 std::span<T> u1_array = u.x()->mutable_array();
1125 std::span<const T> u0_array = v.x()->array();
1126 std::copy(u0_array.begin(), u0_array.end(), u1_array.begin());
1127 }
1128 else
1129 {
1130 // Get mesh and check that functions share the same mesh
1131 if (auto mesh_v = v.function_space()->mesh(); mesh != mesh_v)
1132 {
1133 impl::interpolate_nonmatching_meshes(u, v, cells, nmm_interpolation_data);
1134 return;
1135 }
1136 else
1137 {
1138
1139 // Get elements and check value shape
1140 auto fs0 = v.function_space();
1141 auto element0 = fs0->element();
1142 assert(element0);
1143 auto fs1 = u.function_space();
1144 auto element1 = fs1->element();
1145 assert(element1);
1146 if (fs0->value_shape().size() != fs1->value_shape().size()
1147 or !std::equal(fs0->value_shape().begin(), fs0->value_shape().end(),
1148 fs1->value_shape().begin()))
1149 {
1150 throw std::runtime_error(
1151 "Interpolation: elements have different value dimensions");
1152 }
1153
1154 if (element1 == element0 or *element1 == *element0)
1155 {
1156 // Same element, different dofmaps (or just a subset of cells)
1157
1158 const int tdim = mesh->topology()->dim();
1159 auto cell_map = mesh->topology()->index_map(tdim);
1160 assert(cell_map);
1161
1162 assert(element1->block_size() == element0->block_size());
1163
1164 // Get dofmaps
1165 std::shared_ptr<const DofMap> dofmap0 = v.function_space()->dofmap();
1166 assert(dofmap0);
1167 std::shared_ptr<const DofMap> dofmap1 = u.function_space()->dofmap();
1168 assert(dofmap1);
1169
1170 std::span<T> u1_array = u.x()->mutable_array();
1171 std::span<const T> u0_array = v.x()->array();
1172
1173 // Iterate over mesh and interpolate on each cell
1174 const int bs0 = dofmap0->bs();
1175 const int bs1 = dofmap1->bs();
1176 for (auto c : cells)
1177 {
1178 std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
1179 std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
1180 assert(bs0 * dofs0.size() == bs1 * dofs1.size());
1181 for (std::size_t i = 0; i < dofs0.size(); ++i)
1182 {
1183 for (int k = 0; k < bs0; ++k)
1184 {
1185 int index = bs0 * i + k;
1186 std::div_t dv1 = std::div(index, bs1);
1187 u1_array[bs1 * dofs1[dv1.quot] + dv1.rem]
1188 = u0_array[bs0 * dofs0[i] + k];
1189 }
1190 }
1191 }
1192 }
1193 else if (element1->map_type() == element0->map_type())
1194 {
1195 // Different elements, same basis function map type
1196 impl::interpolate_same_map(u, v, cells);
1197 }
1198 else
1199 {
1200 // Different elements with different maps for basis functions
1201 impl::interpolate_nonmatching_maps(u, v, cells);
1202 }
1203 }
1204 }
1205}
1206} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
A CoordinateElement manages coordinate mappings for isoparametric cells.
Definition CoordinateElement.h:38
static void compute_jacobian(const U &dphi, const V &cell_geometry, W &&J)
Compute Jacobian for a cell with given geometry using the basis functions and first order derivatives...
Definition CoordinateElement.h:105
void tabulate(int nd, std::span< const T > X, std::array< std::size_t, 2 > shape, std::span< T > basis) const
Evaluate basis values and derivatives at set of points.
Definition CoordinateElement.cpp:53
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition CoordinateElement.h:114
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Shape of array to fill when calling tabulate.
Definition CoordinateElement.cpp:46
int dim() const
The dimension of the coordinate element space.
Definition CoordinateElement.cpp:193
static double compute_jacobian_determinant(const U &J, std::span< typename U::value_type > w)
Compute the determinant of the Jacobian.
Definition CoordinateElement.h:131
Finite Element, containing the dof layout on a reference element, and various methods for evaluating ...
Definition FiniteElement.h:29
std::pair< std::vector< geometry_type >, std::array< std::size_t, 2 > > interpolation_points() const
Points on the reference cell at which an expression needs to be evaluated in order to interpolate the...
Definition FiniteElement.cpp:485
This class represents a function in a finite element function space , given by.
Definition Function.h:45
Geometry stores the geometry imposed on a mesh.
Definition Geometry.h:33
const fem::CoordinateElement< value_type > & cmap() const
The element that describes the geometry map.
Definition Geometry.h:180
int dim() const
Return Euclidean dimension of coordinate system.
Definition Geometry.h:116
MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const std::int32_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > dofmap() const
DofMap for the geometry.
Definition Geometry.h:123
std::span< const value_type > x() const
Access geometry degrees-of-freedom data (const version).
Definition Geometry.h:168
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
std::shared_ptr< Topology > topology()
Get mesh topology.
Definition Mesh.h:64
Geometry< T > & geometry()
Get mesh geometry.
Definition Mesh.h:76
Definition interpolate.h:32
void cells(la::SparsityPattern &pattern, std::array< std::span< const std::int32_t >, 2 > cells, std::array< std::reference_wrapper< const DofMap >, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition sparsitybuild.cpp:15
Finite element method functionality.
Definition assemble_matrix_impl.h:26
std::tuple< std::vector< std::int32_t >, std::vector< std::int32_t >, std::vector< T >, std::vector< std::int32_t > > create_nonmatching_meshes_interpolation_data(const mesh::Geometry< T > &geometry0, const FiniteElement< T > &element0, const mesh::Mesh< T > &mesh1, std::span< const std::int32_t > cells, T padding)
Generate data needed to interpolate discrete functions across different meshes.
Definition interpolate.h:1046
void interpolate(Function< T, U > &u, std::span< const T > f, std::array< std::size_t, 2 > fshape, std::span< const std::int32_t > cells)
Interpolate an expression f(x) in a finite element space.
Definition interpolate.h:727
std::vector< T > interpolation_coords(const fem::FiniteElement< T > &element, const mesh::Geometry< T > &geometry, std::span< const std::int32_t > cells)
Compute the evaluation points in the physical space at which an expression should be computed to inte...
Definition interpolate.h:54