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