DOLFINx
DOLFINx C++ interface
interpolate.h
1// Copyright (C) 2020-2021 Garth N. Wells, Igor A. Baratta
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 <dolfinx/common/IndexMap.h>
15#include <dolfinx/mesh/Mesh.h>
16#include <functional>
17#include <numeric>
18#include <span>
19#include <vector>
20
21namespace dolfinx::fem
22{
23template <typename T>
24class Function;
25
26namespace impl
27{
36template <typename U, typename V, typename T>
37void interpolation_apply(const U& Pi, const V& data, std::span<T> coeffs,
38 int bs)
39{
40 static_assert(U::rank() == 2, "Must be rank 2");
41 static_assert(V::rank() == 2, "Must be rank 2");
42
43 // Compute coefficients = Pi * x (matrix-vector multiply)
44 if (bs == 1)
45 {
46 assert(data.extent(0) * data.extent(1) == Pi.extent(1));
47 for (std::size_t i = 0; i < Pi.extent(0); ++i)
48 {
49 coeffs[i] = 0.0;
50 for (std::size_t k = 0; k < data.extent(1); ++k)
51 for (std::size_t j = 0; j < data.extent(0); ++j)
52 coeffs[i] += Pi(i, k * data.extent(0) + j) * data(j, k);
53 }
54 }
55 else
56 {
57 const std::size_t cols = Pi.extent(1);
58 assert(data.extent(0) == Pi.extent(1));
59 assert(data.extent(1) == bs);
60 for (int k = 0; k < bs; ++k)
61 {
62 for (std::size_t i = 0; i < Pi.extent(0); ++i)
63 {
64 T acc = 0;
65 for (std::size_t j = 0; j < cols; ++j)
66 acc += Pi(i, j) * data(j, k);
67 coeffs[bs * i + k] = acc;
68 }
69 }
70 }
71}
72
83template <typename T>
84void interpolate_same_map(Function<T>& u1, const Function<T>& u0,
85 std::span<const std::int32_t> cells)
86{
87 auto V0 = u0.function_space();
88 assert(V0);
89 auto V1 = u1.function_space();
90 assert(V1);
91 auto mesh = V0->mesh();
92 assert(mesh);
93
94 std::shared_ptr<const FiniteElement> element0 = V0->element();
95 assert(element0);
96 std::shared_ptr<const FiniteElement> element1 = V1->element();
97 assert(element1);
98
99 const int tdim = mesh->topology().dim();
100 auto map = mesh->topology().index_map(tdim);
101 assert(map);
102 std::span<T> u1_array = u1.x()->mutable_array();
103 std::span<const T> u0_array = u0.x()->array();
104
105 std::span<const std::uint32_t> cell_info;
106 if (element1->needs_dof_transformations()
107 or element0->needs_dof_transformations())
108 {
109 mesh->topology_mutable().create_entity_permutations();
110 cell_info = std::span(mesh->topology().get_cell_permutation_info());
111 }
112
113 // Get dofmaps
114 auto dofmap1 = V1->dofmap();
115 auto dofmap0 = V0->dofmap();
116
117 // Get block sizes and dof transformation operators
118 const int bs1 = dofmap1->bs();
119 const int bs0 = dofmap0->bs();
120 auto apply_dof_transformation
121 = element0->get_dof_transformation_function<T>(false, true, false);
122 auto apply_inverse_dof_transform
123 = element1->get_dof_transformation_function<T>(true, true, false);
124
125 // Creat working array
126 std::vector<T> local0(element0->space_dimension());
127 std::vector<T> local1(element1->space_dimension());
128
129 // Create interpolation operator
130 const auto [i_m, im_shape]
131 = element1->create_interpolation_operator(*element0);
132
133 // Iterate over mesh and interpolate on each cell
134 for (auto c : cells)
135 {
136 std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
137 for (std::size_t i = 0; i < dofs0.size(); ++i)
138 for (int k = 0; k < bs0; ++k)
139 local0[bs0 * i + k] = u0_array[bs0 * dofs0[i] + k];
140
141 apply_dof_transformation(local0, cell_info, c, 1);
142
143 // FIXME: Get compile-time ranges from Basix
144 // Apply interpolation operator
145 std::fill(local1.begin(), local1.end(), 0);
146 for (std::size_t i = 0; i < im_shape[0]; ++i)
147 for (std::size_t j = 0; j < im_shape[1]; ++j)
148 local1[i] += i_m[im_shape[1] * i + j] * local0[j];
149
150 apply_inverse_dof_transform(local1, cell_info, c, 1);
151
152 std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
153 for (std::size_t i = 0; i < dofs1.size(); ++i)
154 for (int k = 0; k < bs1; ++k)
155 u1_array[bs1 * dofs1[i] + k] = local1[bs1 * i + k];
156 }
157}
158
168template <typename T>
169void interpolate_nonmatching_maps(Function<T>& u1, const Function<T>& u0,
170 std::span<const std::int32_t> cells)
171{
172 // Get mesh
173 auto V0 = u0.function_space();
174 assert(V0);
175 auto mesh = V0->mesh();
176 assert(mesh);
177
178 // Mesh dims
179 const int tdim = mesh->topology().dim();
180 const int gdim = mesh->geometry().dim();
181
182 // Get elements
183 auto V1 = u1.function_space();
184 assert(V1);
185 std::shared_ptr<const FiniteElement> element0 = V0->element();
186 assert(element0);
187 std::shared_ptr<const FiniteElement> element1 = V1->element();
188 assert(element1);
189
190 std::span<const std::uint32_t> cell_info;
191 if (element1->needs_dof_transformations()
192 or element0->needs_dof_transformations())
193 {
194 mesh->topology_mutable().create_entity_permutations();
195 cell_info = std::span(mesh->topology().get_cell_permutation_info());
196 }
197
198 // Get dofmaps
199 auto dofmap0 = V0->dofmap();
200 auto dofmap1 = V1->dofmap();
201
202 const auto [X, Xshape] = element1->interpolation_points();
203
204 // Get block sizes and dof transformation operators
205 const int bs0 = element0->block_size();
206 const int bs1 = element1->block_size();
207 const auto apply_dof_transformation0
208 = element0->get_dof_transformation_function<double>(false, false, false);
209 const auto apply_inverse_dof_transform1
210 = element1->get_dof_transformation_function<T>(true, true, false);
211
212 // Get sizes of elements
213 const std::size_t dim0 = element0->space_dimension() / bs0;
214 const std::size_t value_size_ref0 = element0->reference_value_size() / bs0;
215 const std::size_t value_size0 = element0->value_size() / bs0;
216
217 // Get geometry data
218 const CoordinateElement& cmap = mesh->geometry().cmap();
219 const graph::AdjacencyList<std::int32_t>& x_dofmap
220 = mesh->geometry().dofmap();
221 const std::size_t num_dofs_g = cmap.dim();
222 std::span<const double> x_g = mesh->geometry().x();
223
224 namespace stdex = std::experimental;
225 using cmdspan2_t
226 = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
227 using cmdspan4_t
228 = stdex::mdspan<const double, stdex::dextents<std::size_t, 4>>;
229 using mdspan2_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
230 using mdspan3_t = stdex::mdspan<double, stdex::dextents<std::size_t, 3>>;
231 using mdspan3T_t = stdex::mdspan<T, stdex::dextents<std::size_t, 3>>;
232
233 // Evaluate coordinate map basis at reference interpolation points
234 const std::array<std::size_t, 4> phi_shape
235 = cmap.tabulate_shape(1, Xshape[0]);
236 std::vector<double> phi_b(
237 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
238 cmdspan4_t phi(phi_b.data(), phi_shape);
239 cmap.tabulate(1, X, Xshape, phi_b);
240
241 // Evaluate v basis functions at reference interpolation points
242 const auto [_basis_derivatives_reference0, b0shape]
243 = element0->tabulate(X, Xshape, 0);
244 cmdspan4_t basis_derivatives_reference0(_basis_derivatives_reference0.data(),
245 b0shape);
246
247 // Create working arrays
248 std::vector<T> local1(element1->space_dimension());
249 std::vector<T> coeffs0(element0->space_dimension());
250
251 std::vector<double> basis0_b(Xshape[0] * dim0 * value_size0);
252 mdspan3_t basis0(basis0_b.data(), Xshape[0], dim0, value_size0);
253
254 std::vector<double> basis_reference0_b(Xshape[0] * dim0 * value_size_ref0);
255 mdspan3_t basis_reference0(basis_reference0_b.data(), Xshape[0], dim0,
256 value_size_ref0);
257
258 std::vector<T> values0_b(Xshape[0] * 1 * element1->value_size());
259 mdspan3T_t values0(values0_b.data(), Xshape[0], 1, element1->value_size());
260
261 std::vector<T> mapped_values_b(Xshape[0] * 1 * element1->value_size());
262 mdspan3T_t mapped_values0(mapped_values_b.data(), Xshape[0], 1,
263 element1->value_size());
264
265 std::vector<double> coord_dofs_b(num_dofs_g * gdim);
266 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
267
268 std::vector<double> J_b(Xshape[0] * gdim * tdim);
269 mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
270 std::vector<double> K_b(Xshape[0] * tdim * gdim);
271 mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
272 std::vector<double> detJ(Xshape[0]);
273 std::vector<double> det_scratch(2 * gdim * tdim);
274
275 // Get interpolation operator
276 const auto [_Pi_1, pi_shape] = element1->interpolation_operator();
277 cmdspan2_t Pi_1(_Pi_1.data(), pi_shape);
278
279 using u_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
280 using U_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
281 using J_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
282 using K_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
283 auto push_forward_fn0
284 = element0->basix_element().map_fn<u_t, U_t, J_t, K_t>();
285
286 using v_t = stdex::mdspan<const T, stdex::dextents<std::size_t, 2>>;
287 using V_t = stdex::mdspan<T, stdex::dextents<std::size_t, 2>>;
288 auto pull_back_fn1 = element1->basix_element().map_fn<V_t, v_t, K_t, J_t>();
289
290 // Iterate over mesh and interpolate on each cell
291 std::span<const T> array0 = u0.x()->array();
292 std::span<T> array1 = u1.x()->mutable_array();
293 for (auto c : cells)
294 {
295 // Get cell geometry (coordinate dofs)
296 auto x_dofs = x_dofmap.links(c);
297 for (std::size_t i = 0; i < num_dofs_g; ++i)
298 {
299 const int pos = 3 * x_dofs[i];
300 for (std::size_t j = 0; j < gdim; ++j)
301 coord_dofs(i, j) = x_g[pos + j];
302 }
303
304 // Compute Jacobians and reference points for current cell
305 std::fill(J_b.begin(), J_b.end(), 0);
306 for (std::size_t p = 0; p < Xshape[0]; ++p)
307 {
308 auto dphi = stdex::submdspan(phi, std::pair(1, tdim + 1), p,
309 stdex::full_extent, 0);
310
311 auto _J = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
312 cmap.compute_jacobian(dphi, coord_dofs, _J);
313 auto _K = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
314 cmap.compute_jacobian_inverse(_J, _K);
315 detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
316 }
317
318 // Copy evaluated basis on reference, apply DOF transformations, and
319 // push forward to physical element
320 for (std::size_t k0 = 0; k0 < basis_reference0.extent(0); ++k0)
321 for (std::size_t k1 = 0; k1 < basis_reference0.extent(1); ++k1)
322 for (std::size_t k2 = 0; k2 < basis_reference0.extent(2); ++k2)
323 basis_reference0(k0, k1, k2)
324 = basis_derivatives_reference0(0, k0, k1, k2);
325
326 for (std::size_t p = 0; p < Xshape[0]; ++p)
327 {
328 apply_dof_transformation0(
329 std::span(basis_reference0_b.data() + p * dim0 * value_size_ref0,
330 dim0 * value_size_ref0),
331 cell_info, c, value_size_ref0);
332 }
333
334 for (std::size_t i = 0; i < basis0.extent(0); ++i)
335 {
336 auto _u
337 = stdex::submdspan(basis0, i, stdex::full_extent, stdex::full_extent);
338 auto _U = stdex::submdspan(basis_reference0, i, stdex::full_extent,
339 stdex::full_extent);
340 auto _K = stdex::submdspan(K, i, stdex::full_extent, stdex::full_extent);
341 auto _J = stdex::submdspan(J, i, stdex::full_extent, stdex::full_extent);
342 push_forward_fn0(_u, _U, _J, detJ[i], _K);
343 }
344
345 // Copy expansion coefficients for v into local array
346 const int dof_bs0 = dofmap0->bs();
347 std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
348 for (std::size_t i = 0; i < dofs0.size(); ++i)
349 for (int k = 0; k < dof_bs0; ++k)
350 coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k];
351
352 // Evaluate v at the interpolation points (physical space values)
353 for (std::size_t p = 0; p < Xshape[0]; ++p)
354 {
355 for (int k = 0; k < bs0; ++k)
356 {
357 for (std::size_t j = 0; j < value_size0; ++j)
358 {
359 T acc = 0;
360 for (std::size_t i = 0; i < dim0; ++i)
361 acc += coeffs0[bs0 * i + k] * basis0(p, i, j);
362 values0(p, 0, j * bs0 + k) = acc;
363 }
364 }
365 }
366
367 // Pull back the physical values to the u reference
368 for (std::size_t i = 0; i < values0.extent(0); ++i)
369 {
370 auto _u = stdex::submdspan(values0, i, stdex::full_extent,
371 stdex::full_extent);
372 auto _U = stdex::submdspan(mapped_values0, i, stdex::full_extent,
373 stdex::full_extent);
374 auto _K = stdex::submdspan(K, i, stdex::full_extent, stdex::full_extent);
375 auto _J = stdex::submdspan(J, i, stdex::full_extent, stdex::full_extent);
376 pull_back_fn1(_U, _u, _K, 1.0 / detJ[i], _J);
377 }
378
379 auto values = stdex::submdspan(mapped_values0, stdex::full_extent, 0,
380 stdex::full_extent);
381 interpolation_apply(Pi_1, values, std::span(local1), bs1);
382 apply_inverse_dof_transform1(local1, cell_info, c, 1);
383
384 // Copy local coefficients to the correct position in u dof array
385 const int dof_bs1 = dofmap1->bs();
386 std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
387 for (std::size_t i = 0; i < dofs1.size(); ++i)
388 for (int k = 0; k < dof_bs1; ++k)
389 array1[dof_bs1 * dofs1[i] + k] = local1[dof_bs1 * i + k];
390 }
391}
392} // namespace impl
393
404std::vector<double> interpolation_coords(const FiniteElement& element,
405 const mesh::Mesh& mesh,
406 std::span<const std::int32_t> cells);
407
420template <typename T>
421void interpolate(Function<T>& u, std::span<const T> f,
422 std::array<std::size_t, 2> fshape,
423 std::span<const std::int32_t> cells)
424{
425 namespace stdex = std::experimental;
426 using cmdspan2_t
427 = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
428 using cmdspan4_t
429 = stdex::mdspan<const double, stdex::dextents<std::size_t, 4>>;
430 using mdspan2_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
431 using mdspan3_t = stdex::mdspan<double, stdex::dextents<std::size_t, 3>>;
432
433 const std::shared_ptr<const FiniteElement> element
434 = u.function_space()->element();
435 assert(element);
436 const int element_bs = element->block_size();
437 if (int num_sub = element->num_sub_elements();
438 num_sub > 0 and num_sub != element_bs)
439 {
440 throw std::runtime_error("Cannot directly interpolate a mixed space. "
441 "Interpolate into subspaces.");
442 }
443
444 if (fshape[0] != (std::size_t)element->value_size())
445 throw std::runtime_error("Interpolation data has the wrong shape/size.");
446
447 // Get mesh
448 assert(u.function_space());
449 auto mesh = u.function_space()->mesh();
450 assert(mesh);
451
452 const int gdim = mesh->geometry().dim();
453 const int tdim = mesh->topology().dim();
454
455 std::span<const std::uint32_t> cell_info;
456 if (element->needs_dof_transformations())
457 {
458 mesh->topology_mutable().create_entity_permutations();
459 cell_info = std::span(mesh->topology().get_cell_permutation_info());
460 }
461
462 const std::size_t f_shape1 = f.size() / element->value_size();
463 stdex::mdspan<const T, stdex::dextents<std::size_t, 2>> _f(f.data(), fshape);
464
465 // Get dofmap
466 const auto dofmap = u.function_space()->dofmap();
467 assert(dofmap);
468 const int dofmap_bs = dofmap->bs();
469
470 // Loop over cells and compute interpolation dofs
471 const int num_scalar_dofs = element->space_dimension() / element_bs;
472 const int value_size = element->value_size() / element_bs;
473
474 std::span<T> coeffs = u.x()->mutable_array();
475 std::vector<T> _coeffs(num_scalar_dofs);
476
477 // This assumes that any element with an identity interpolation matrix
478 // is a point evaluation
479 if (element->map_ident() && element->interpolation_ident())
480 {
481 // Point evaluation element *and* the geometric map is the identity,
482 // e.g. not Piola mapped
483
484 auto apply_inv_transpose_dof_transformation
485 = element->get_dof_transformation_function<T>(true, true, true);
486
487 // Loop over cells
488 for (std::size_t c = 0; c < cells.size(); ++c)
489 {
490 const std::int32_t cell = cells[c];
491 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
492 for (int k = 0; k < element_bs; ++k)
493 {
494 // num_scalar_dofs is the number of interpolation points per
495 // cell in this case (interpolation matrix is identity)
496 std::copy_n(std::next(f.begin(), k * f_shape1 + c * num_scalar_dofs),
497 num_scalar_dofs, _coeffs.begin());
498 apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
499 for (int i = 0; i < num_scalar_dofs; ++i)
500 {
501 const int dof = i * element_bs + k;
502 std::div_t pos = std::div(dof, dofmap_bs);
503 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
504 }
505 }
506 }
507 }
508 else if (element->map_ident())
509 {
510 // Not a point evaluation, but the geometric map is the identity,
511 // e.g. not Piola mapped
512
513 const int element_vs = element->value_size() / element_bs;
514
515 if (element_vs > 1 && element_bs > 1)
516 {
517 throw std::runtime_error(
518 "Interpolation into this element not supported.");
519 }
520
521 // Get interpolation operator
522 const auto [_Pi, pi_shape] = element->interpolation_operator();
523 cmdspan2_t Pi(_Pi.data(), pi_shape);
524 const std::size_t num_interp_points = Pi.extent(1);
525 assert(Pi.extent(0) == num_scalar_dofs);
526
527 auto apply_inv_transpose_dof_transformation
528 = element->get_dof_transformation_function<T>(true, true, true);
529
530 // Loop over cells
531 std::vector<T> ref_data_b(num_interp_points);
532 stdex::mdspan<T, stdex::extents<std::size_t, stdex::dynamic_extent, 1>>
533 ref_data(ref_data_b.data(), num_interp_points, 1);
534 for (std::size_t c = 0; c < cells.size(); ++c)
535 {
536 const std::int32_t cell = cells[c];
537 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
538 for (int k = 0; k < element_bs; ++k)
539 {
540 for (int i = 0; i < element_vs; ++i)
541 {
542 std::copy_n(
543 std::next(f.begin(), (i + k) * f_shape1
544 + c * num_interp_points / element_vs),
545 num_interp_points / element_vs,
546 std::next(ref_data_b.begin(),
547 i * num_interp_points / element_vs));
548 }
549 impl::interpolation_apply(Pi, ref_data, std::span(_coeffs), 1);
550 apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
551 for (int i = 0; i < num_scalar_dofs; ++i)
552 {
553 const int dof = i * element_bs + k;
554 std::div_t pos = std::div(dof, dofmap_bs);
555 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
556 }
557 }
558 }
559 }
560 else
561 {
562 // Get the interpolation points on the reference cells
563 const auto [X, Xshape] = element->interpolation_points();
564 if (X.empty())
565 {
566 throw std::runtime_error(
567 "Interpolation into this space is not yet supported.");
568 }
569
570 if (_f.extent(1) != cells.size() * Xshape[0])
571 throw std::runtime_error("Interpolation data has the wrong shape.");
572
573 // Get coordinate map
574 const CoordinateElement& cmap = mesh->geometry().cmap();
575
576 // Get geometry data
578 = mesh->geometry().dofmap();
579 const int num_dofs_g = cmap.dim();
580 std::span<const double> x_g = mesh->geometry().x();
581
582 // Create data structures for Jacobian info
583 std::vector<double> J_b(Xshape[0] * gdim * tdim);
584 mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
585 std::vector<double> K_b(Xshape[0] * tdim * gdim);
586 mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
587 std::vector<double> detJ(Xshape[0]);
588 std::vector<double> det_scratch(2 * gdim * tdim);
589
590 std::vector<double> coord_dofs_b(num_dofs_g * gdim);
591 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
592
593 std::vector<T> ref_data_b(Xshape[0] * 1 * value_size);
594 stdex::mdspan<T, stdex::dextents<std::size_t, 3>> ref_data(
595 ref_data_b.data(), Xshape[0], 1, value_size);
596
597 std::vector<T> _vals_b(Xshape[0] * 1 * value_size);
598 stdex::mdspan<T, stdex::dextents<std::size_t, 3>> _vals(
599 _vals_b.data(), Xshape[0], 1, value_size);
600
601 // Tabulate 1st derivative of shape functions at interpolation
602 // coords
603 std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(1, Xshape[0]);
604 std::vector<double> phi_b(
605 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
606 cmdspan4_t phi(phi_b.data(), phi_shape);
607 cmap.tabulate(1, X, Xshape, phi_b);
608 auto dphi = stdex::submdspan(phi, std::pair(1, tdim + 1),
609 stdex::full_extent, stdex::full_extent, 0);
610
611 const std::function<void(const std::span<T>&,
612 const std::span<const std::uint32_t>&,
613 std::int32_t, int)>
614 apply_inverse_transpose_dof_transformation
615 = element->get_dof_transformation_function<T>(true, true);
616
617 // Get interpolation operator
618 const auto [_Pi, pi_shape] = element->interpolation_operator();
619 cmdspan2_t Pi(_Pi.data(), pi_shape);
620
621 namespace stdex = std::experimental;
622 using u_t = stdex::mdspan<const T, stdex::dextents<std::size_t, 2>>;
623 using U_t = stdex::mdspan<T, stdex::dextents<std::size_t, 2>>;
624 using J_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
625 using K_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
626 auto pull_back_fn = element->basix_element().map_fn<U_t, u_t, J_t, K_t>();
627
628 for (std::size_t c = 0; c < cells.size(); ++c)
629 {
630 const std::int32_t cell = cells[c];
631 auto x_dofs = x_dofmap.links(cell);
632 for (int i = 0; i < num_dofs_g; ++i)
633 {
634 const int pos = 3 * x_dofs[i];
635 for (int j = 0; j < gdim; ++j)
636 coord_dofs(i, j) = x_g[pos + j];
637 }
638
639 // Compute J, detJ and K
640 std::fill(J_b.begin(), J_b.end(), 0);
641 for (std::size_t p = 0; p < Xshape[0]; ++p)
642 {
643 auto _dphi
644 = stdex::submdspan(dphi, stdex::full_extent, p, stdex::full_extent);
645 auto _J
646 = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
647 cmap.compute_jacobian(_dphi, coord_dofs, _J);
648 auto _K
649 = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
650 cmap.compute_jacobian_inverse(_J, _K);
651 detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
652 }
653
654 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
655 for (int k = 0; k < element_bs; ++k)
656 {
657 // Extract computed expression values for element block k
658 for (int m = 0; m < value_size; ++m)
659 {
660 for (std::size_t k0 = 0; k0 < Xshape[0]; ++k0)
661 {
662 _vals(k0, 0, m)
663 = f[f_shape1 * (k * value_size + m) + c * Xshape[0] + k0];
664 }
665 }
666
667 // Get element degrees of freedom for block
668 for (std::size_t i = 0; i < Xshape[0]; ++i)
669 {
670 auto _u = stdex::submdspan(_vals, i, stdex::full_extent,
671 stdex::full_extent);
672 auto _U = stdex::submdspan(ref_data, i, stdex::full_extent,
673 stdex::full_extent);
674 auto _K
675 = stdex::submdspan(K, i, stdex::full_extent, stdex::full_extent);
676 auto _J
677 = stdex::submdspan(J, i, stdex::full_extent, stdex::full_extent);
678 pull_back_fn(_U, _u, _K, 1.0 / detJ[i], _J);
679 }
680
681 auto ref = stdex::submdspan(ref_data, stdex::full_extent, 0,
682 stdex::full_extent);
683 impl::interpolation_apply(Pi, ref, std::span(_coeffs), element_bs);
684 apply_inverse_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
685
686 // Copy interpolation dofs into coefficient vector
687 assert(_coeffs.size() == num_scalar_dofs);
688 for (int i = 0; i < num_scalar_dofs; ++i)
689 {
690 const int dof = i * element_bs + k;
691 std::div_t pos = std::div(dof, dofmap_bs);
692 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
693 }
694 }
695 }
696 }
697}
698
704template <typename T>
706 const std::span<const std::int32_t>& cells)
707{
708 assert(u.function_space());
709 assert(v.function_space());
710 std::shared_ptr<const mesh::Mesh> mesh = u.function_space()->mesh();
711 assert(mesh);
712
713 auto cell_map0 = mesh->topology().index_map(mesh->topology().dim());
714 assert(cell_map0);
715 std::size_t num_cells0 = cell_map0->size_local() + cell_map0->num_ghosts();
716 if (u.function_space() == v.function_space() and cells.size() == num_cells0)
717 {
718 // Same function spaces and on whole mesh
719 std::span<T> u1_array = u.x()->mutable_array();
720 std::span<const T> u0_array = v.x()->array();
721 std::copy(u0_array.begin(), u0_array.end(), u1_array.begin());
722 }
723 else
724 {
725 // Get mesh and check that functions share the same mesh
726 if (mesh != v.function_space()->mesh())
727 {
728 throw std::runtime_error(
729 "Interpolation on different meshes not supported (yet).");
730 }
731
732 // Get elements and check value shape
733 auto element0 = v.function_space()->element();
734 assert(element0);
735 auto element1 = u.function_space()->element();
736 assert(element1);
737 if (element0->value_shape().size() != element1->value_shape().size()
738 or !std::equal(element0->value_shape().begin(),
739 element0->value_shape().end(),
740 element1->value_shape().begin()))
741 {
742 throw std::runtime_error(
743 "Interpolation: elements have different value dimensions");
744 }
745
746 if (*element1 == *element0)
747 {
748 // Same element, different dofmaps (or just a subset of cells)
749
750 const int tdim = mesh->topology().dim();
751 auto cell_map = mesh->topology().index_map(tdim);
752 assert(cell_map);
753
754 assert(element1->block_size() == element0->block_size());
755
756 // Get dofmaps
757 std::shared_ptr<const DofMap> dofmap0 = v.function_space()->dofmap();
758 assert(dofmap0);
759 std::shared_ptr<const DofMap> dofmap1 = u.function_space()->dofmap();
760 assert(dofmap1);
761
762 std::span<T> u1_array = u.x()->mutable_array();
763 std::span<const T> u0_array = v.x()->array();
764
765 // Iterate over mesh and interpolate on each cell
766 const int bs0 = dofmap0->bs();
767 const int bs1 = dofmap1->bs();
768 for (auto c : cells)
769 {
770 std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
771 std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
772 assert(bs0 * dofs0.size() == bs1 * dofs1.size());
773 for (std::size_t i = 0; i < dofs0.size(); ++i)
774 {
775 for (int k = 0; k < bs0; ++k)
776 {
777 int index = bs0 * i + k;
778 std::div_t dv1 = std::div(index, bs1);
779 u1_array[bs1 * dofs1[dv1.quot] + dv1.rem]
780 = u0_array[bs0 * dofs0[i] + k];
781 }
782 }
783 }
784 }
785 else if (element1->map_type() == element0->map_type())
786 {
787 // Different elements, same basis function map type
788 impl::interpolate_same_map(u, v, cells);
789 }
790 else
791 {
792 // Different elements with different maps for basis functions
793 impl::interpolate_nonmatching_maps(u, v, cells);
794 }
795 }
796}
797
798} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
A CoordinateElement manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:32
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:100
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition: CoordinateElement.h:109
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 FiniteElement::tabulate
Definition: CoordinateElement.cpp:42
void tabulate(int nd, std::span< const double > X, std::array< std::size_t, 2 > shape, std::span< double > basis) const
Evaluate basis values and derivatives at set of points.
Definition: CoordinateElement.cpp:47
int dim() const
The dimension of the geometry element space.
Definition: CoordinateElement.cpp:183
static double compute_jacobian_determinant(const U &J, std::span< typename U::value_type > w)
Compute the determinant of the Jacobian.
Definition: CoordinateElement.h:126
This class represents a function in a finite element function space , given by.
Definition: Function.h:45
std::shared_ptr< const FunctionSpace > function_space() const
Access the function space.
Definition: Function.h:136
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:142
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:26
std::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:111
Finite element method functionality.
Definition: assemble_matrix_impl.h:25
void interpolate(Function< T > &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:421
std::vector< double > interpolation_coords(const FiniteElement &element, const mesh::Mesh &mesh, 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.cpp:16