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 <vector>
19#include <xtensor/xarray.hpp>
20#include <xtensor/xtensor.hpp>
21#include <xtensor/xview.hpp>
22#include <xtl/xspan.hpp>
23
24namespace dolfinx::fem
25{
26template <typename T>
27class Function;
28
29namespace impl
30{
39template <typename U, typename V, typename T>
40void interpolation_apply(const U& Pi, const V& data, std::vector<T>& coeffs,
41 int bs)
42{
43 // Compute coefficients = Pi * x (matrix-vector multiply)
44 if (bs == 1)
45 {
46 assert(data.shape(0) * data.shape(1) == Pi.shape(1));
47 for (std::size_t i = 0; i < Pi.shape(0); ++i)
48 {
49 coeffs[i] = 0.0;
50 for (std::size_t k = 0; k < data.shape(1); ++k)
51 for (std::size_t j = 0; j < data.shape(0); ++j)
52 coeffs[i] += Pi(i, k * data.shape(0) + j) * data(j, k);
53 }
54 }
55 else
56 {
57 const std::size_t cols = Pi.shape(1);
58 assert(data.shape(0) == Pi.shape(1));
59 assert(data.shape(1) == bs);
60 for (int k = 0; k < bs; ++k)
61 {
62 for (std::size_t i = 0; i < Pi.shape(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 const xtl::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 xtl::span<T> u1_array = u1.x()->mutable_array();
103 xtl::span<const T> u0_array = u0.x()->array();
104
105 xtl::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 = xtl::span(mesh->topology().get_cell_permutation_info());
111 }
112
113 // Get dofmaps
114 auto dofmap1 = V1->dofmap();
115 auto dofmap0 = V0->dofmap();
116
117 // Create interpolation operator
118 const xt::xtensor<double, 2> i_m
119 = element1->create_interpolation_operator(*element0);
120
121 // Get block sizes and dof transformation operators
122 const int bs1 = dofmap1->bs();
123 const int bs0 = dofmap0->bs();
124 auto apply_dof_transformation
125 = element0->get_dof_transformation_function<T>(false, true, false);
126 auto apply_inverse_dof_transform
127 = element1->get_dof_transformation_function<T>(true, true, false);
128
129 // Creat working array
130 std::vector<T> local0(element0->space_dimension());
131 std::vector<T> local1(element1->space_dimension());
132
133 // Iterate over mesh and interpolate on each cell
134 for (auto c : cells)
135 {
136 xtl::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 < i_m.shape(0); ++i)
147 for (std::size_t j = 0; j < i_m.shape(1); ++j)
148 local1[i] += i_m(i, j) * local0[j];
149
150 apply_inverse_dof_transform(local1, cell_info, c, 1);
151
152 xtl::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 const xtl::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 xtl::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 = xtl::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 xt::xtensor<double, 2> X = 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 fem::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 xtl::span<const double> x_g = mesh->geometry().x();
223
224 // Evaluate coordinate map basis at reference interpolation points
225 xt::xtensor<double, 4> phi(cmap.tabulate_shape(1, X.shape(0)));
226 xt::xtensor<double, 2> dphi;
227 cmap.tabulate(1, X, phi);
228 dphi = xt::view(phi, xt::range(1, tdim + 1), 0, xt::all(), 0);
229
230 // Evaluate v basis functions at reference interpolation points
231 xt::xtensor<double, 4> basis_derivatives_reference0(
232 {1, X.shape(0), dim0, value_size_ref0});
233 element0->tabulate(basis_derivatives_reference0, X, 0);
234
235 // Create working arrays
236 std::vector<T> local1(element1->space_dimension());
237 std::vector<T> coeffs0(element0->space_dimension());
238 xt::xtensor<double, 3> basis0({X.shape(0), dim0, value_size0});
239 xt::xtensor<double, 3> basis_reference0({X.shape(0), dim0, value_size_ref0});
240 xt::xtensor<T, 3> values0({X.shape(0), 1, element1->value_size()});
241 xt::xtensor<T, 3> mapped_values0({X.shape(0), 1, element1->value_size()});
242 xt::xtensor<double, 2> coordinate_dofs({num_dofs_g, gdim});
243 xt::xtensor<double, 3> J({X.shape(0), gdim, tdim});
244 xt::xtensor<double, 3> K({X.shape(0), tdim, gdim});
245 std::vector<double> detJ(X.shape(0));
246
247 // Get interpolation operator
248 const xt::xtensor<double, 2>& Pi_1 = element1->interpolation_operator();
249
250 using u_t = xt::xview<decltype(basis_reference0)&, std::size_t,
251 xt::xall<std::size_t>, xt::xall<std::size_t>>;
252 using U_t = xt::xview<decltype(basis_reference0)&, std::size_t,
253 xt::xall<std::size_t>, xt::xall<std::size_t>>;
254 using J_t = xt::xview<decltype(J)&, std::size_t, xt::xall<std::size_t>,
255 xt::xall<std::size_t>>;
256 using K_t = xt::xview<decltype(K)&, std::size_t, xt::xall<std::size_t>,
257 xt::xall<std::size_t>>;
258 auto push_forward_fn0 = element0->map_fn<u_t, U_t, J_t, K_t>();
259
260 using u1_t = xt::xview<decltype(values0)&, std::size_t, xt::xall<std::size_t>,
261 xt::xall<std::size_t>>;
262 using U1_t = xt::xview<decltype(mapped_values0)&, std::size_t,
263 xt::xall<std::size_t>, xt::xall<std::size_t>>;
264 auto pull_back_fn1 = element1->map_fn<U1_t, u1_t, K_t, J_t>();
265
266 // Iterate over mesh and interpolate on each cell
267 xtl::span<const T> array0 = u0.x()->array();
268 xtl::span<T> array1 = u1.x()->mutable_array();
269 for (auto c : cells)
270 {
271 // Get cell geometry (coordinate dofs)
272 auto x_dofs = x_dofmap.links(c);
273 for (std::size_t i = 0; i < num_dofs_g; ++i)
274 {
275 const int pos = 3 * x_dofs[i];
276 for (std::size_t j = 0; j < gdim; ++j)
277 coordinate_dofs(i, j) = x_g[pos + j];
278 }
279
280 // Compute Jacobians and reference points for current cell
281 J.fill(0);
282 for (std::size_t p = 0; p < X.shape(0); ++p)
283 {
284 auto _J = xt::view(J, p, xt::all(), xt::all());
285 cmap.compute_jacobian(dphi, coordinate_dofs, _J);
286 cmap.compute_jacobian_inverse(_J, xt::view(K, p, xt::all(), xt::all()));
287 detJ[p] = cmap.compute_jacobian_determinant(_J);
288 }
289
290 // Get evaluated basis on reference, apply DOF transformations, and
291 // push forward to physical element
292 basis_reference0 = xt::view(basis_derivatives_reference0, 0, xt::all(),
293 xt::all(), xt::all());
294 for (std::size_t p = 0; p < X.shape(0); ++p)
295 {
296 apply_dof_transformation0(
297 xtl::span(basis_reference0.data() + p * dim0 * value_size_ref0,
298 dim0 * value_size_ref0),
299 cell_info, c, value_size_ref0);
300 }
301
302 for (std::size_t i = 0; i < basis0.shape(0); ++i)
303 {
304 auto _K = xt::view(K, i, xt::all(), xt::all());
305 auto _J = xt::view(J, i, xt::all(), xt::all());
306 auto _u = xt::view(basis0, i, xt::all(), xt::all());
307 auto _U = xt::view(basis_reference0, i, xt::all(), xt::all());
308 push_forward_fn0(_u, _U, _J, detJ[i], _K);
309 }
310
311 // Copy expansion coefficients for v into local array
312 const int dof_bs0 = dofmap0->bs();
313 xtl::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
314 for (std::size_t i = 0; i < dofs0.size(); ++i)
315 for (int k = 0; k < dof_bs0; ++k)
316 coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k];
317
318 // Evaluate v at the interpolation points (physical space values)
319 for (std::size_t p = 0; p < X.shape(0); ++p)
320 {
321 for (int k = 0; k < bs0; ++k)
322 {
323 for (std::size_t j = 0; j < value_size0; ++j)
324 {
325 T acc = 0;
326 for (std::size_t i = 0; i < dim0; ++i)
327 acc += coeffs0[bs0 * i + k] * basis0(p, i, j);
328 values0(p, 0, j * bs0 + k) = acc;
329 }
330 }
331 }
332
333 // Pull back the physical values to the u reference
334 for (std::size_t i = 0; i < values0.shape(0); ++i)
335 {
336 auto _K = xt::view(K, i, xt::all(), xt::all());
337 auto _J = xt::view(J, i, xt::all(), xt::all());
338 auto _u = xt::view(values0, i, xt::all(), xt::all());
339 auto _U = xt::view(mapped_values0, i, xt::all(), xt::all());
340 pull_back_fn1(_U, _u, _K, 1.0 / detJ[i], _J);
341 }
342
343 auto _mapped_values0 = xt::view(mapped_values0, xt::all(), 0, xt::all());
344 interpolation_apply(Pi_1, _mapped_values0, local1, bs1);
345 apply_inverse_dof_transform1(local1, cell_info, c, 1);
346
347 // Copy local coefficients to the correct position in u dof array
348 const int dof_bs1 = dofmap1->bs();
349 xtl::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
350 for (std::size_t i = 0; i < dofs1.size(); ++i)
351 for (int k = 0; k < dof_bs1; ++k)
352 array1[dof_bs1 * dofs1[i] + k] = local1[dof_bs1 * i + k];
353 }
354}
355} // namespace impl
356
367xt::xtensor<double, 2>
368interpolation_coords(const fem::FiniteElement& element, const mesh::Mesh& mesh,
369 const xtl::span<const std::int32_t>& cells);
370
382template <typename T>
383void interpolate(Function<T>& u, const xt::xarray<T>& f,
384 const xtl::span<const std::int32_t>& cells)
385{
386 const std::shared_ptr<const FiniteElement> element
387 = u.function_space()->element();
388 assert(element);
389 const int element_bs = element->block_size();
390 if (int num_sub = element->num_sub_elements();
391 num_sub > 0 and num_sub != element_bs)
392 {
393 throw std::runtime_error("Cannot directly interpolate a mixed space. "
394 "Interpolate into subspaces.");
395 }
396
397 // Get mesh
398 assert(u.function_space());
399 auto mesh = u.function_space()->mesh();
400 assert(mesh);
401
402 const int gdim = mesh->geometry().dim();
403 const int tdim = mesh->topology().dim();
404
405 xtl::span<const std::uint32_t> cell_info;
406 if (element->needs_dof_transformations())
407 {
408 mesh->topology_mutable().create_entity_permutations();
409 cell_info = xtl::span(mesh->topology().get_cell_permutation_info());
410 }
411
412 if (f.dimension() == 1)
413 {
414 if (element->value_size() != 1)
415 throw std::runtime_error("Interpolation data has the wrong shape/size.");
416 }
417 else if (f.dimension() == 2)
418 {
419 if (f.shape(0) != element->value_size())
420 throw std::runtime_error("Interpolation data has the wrong shape/size.");
421 }
422 else
423 throw std::runtime_error("Interpolation data has wrong shape.");
424
425 const xtl::span<const T> _f(f.data(), f.size());
426 const std::size_t f_shape1 = _f.size() / element->value_size();
427
428 // Get dofmap
429 const auto dofmap = u.function_space()->dofmap();
430 assert(dofmap);
431 const int dofmap_bs = dofmap->bs();
432
433 // Loop over cells and compute interpolation dofs
434 const int num_scalar_dofs = element->space_dimension() / element_bs;
435 const int value_size = element->value_size() / element_bs;
436
437 xtl::span<T> coeffs = u.x()->mutable_array();
438 std::vector<T> _coeffs(num_scalar_dofs);
439
440 // This assumes that any element with an identity interpolation matrix
441 // is a point evaluation
442 if (element->interpolation_ident())
443 {
444 if (!element->map_ident())
445 throw std::runtime_error("Element does not have identity map.");
446
447 auto apply_inv_transpose_dof_transformation
448 = element->get_dof_transformation_function<T>(true, true, true);
449
450 // Loop over cells
451 for (std::size_t c = 0; c < cells.size(); ++c)
452 {
453 const std::int32_t cell = cells[c];
454 xtl::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
455 for (int k = 0; k < element_bs; ++k)
456 {
457 // num_scalar_dofs is the number of interpolation points per
458 // cell in this case (interpolation matrix is identity)
459 std::copy_n(std::next(_f.begin(), k * f_shape1 + c * num_scalar_dofs),
460 num_scalar_dofs, _coeffs.begin());
461 apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
462 for (int i = 0; i < num_scalar_dofs; ++i)
463 {
464 const int dof = i * element_bs + k;
465 std::div_t pos = std::div(dof, dofmap_bs);
466 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
467 }
468 }
469 }
470 }
471 else if (element->map_ident())
472 {
473 if (f.dimension() != 1)
474 throw std::runtime_error("Interpolation data has the wrong shape.");
475
476 // Get interpolation operator
477 const xt::xtensor<double, 2>& Pi = element->interpolation_operator();
478 const std::size_t num_interp_points = Pi.shape(1);
479 assert(Pi.shape(0) == num_scalar_dofs);
480
481 auto apply_inv_transpose_dof_transformation
482 = element->get_dof_transformation_function<T>(true, true, true);
483
484 // Loop over cells
485 xt::xtensor<T, 2> reference_data({num_interp_points, 1});
486 for (std::size_t c = 0; c < cells.size(); ++c)
487 {
488 const std::int32_t cell = cells[c];
489 xtl::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
490 for (int k = 0; k < element_bs; ++k)
491 {
492 std::copy_n(std::next(_f.begin(), k * f_shape1 + c * num_interp_points),
493 num_interp_points, reference_data.begin());
494
495 impl::interpolation_apply(Pi, reference_data, _coeffs, element_bs);
496
497 apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
498 for (int i = 0; i < num_scalar_dofs; ++i)
499 {
500 const int dof = i * element_bs + k;
501 std::div_t pos = std::div(dof, dofmap_bs);
502 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
503 }
504 }
505 }
506 }
507 else
508 {
509 // Get the interpolation points on the reference cells
510 const xt::xtensor<double, 2>& X = element->interpolation_points();
511 if (X.shape(0) == 0)
512 {
513 throw std::runtime_error(
514 "Interpolation into this space is not yet supported.");
515 }
516
517 if (f.shape(1) != cells.size() * X.shape(0))
518 throw std::runtime_error("Interpolation data has the wrong shape.");
519
520 // Get coordinate map
521 const fem::CoordinateElement& cmap = mesh->geometry().cmap();
522
523 // Get geometry data
525 = mesh->geometry().dofmap();
526 const int num_dofs_g = cmap.dim();
527 xtl::span<const double> x_g = mesh->geometry().x();
528
529 // Create data structures for Jacobian info
530 xt::xtensor<double, 3> J = xt::empty<double>({int(X.shape(0)), gdim, tdim});
531 xt::xtensor<double, 3> K = xt::empty<double>({int(X.shape(0)), tdim, gdim});
532 xt::xtensor<double, 1> detJ = xt::empty<double>({X.shape(0)});
533
534 xt::xtensor<double, 2> coordinate_dofs
535 = xt::empty<double>({num_dofs_g, gdim});
536
537 xt::xtensor<T, 3> reference_data({X.shape(0), 1, value_size});
538 xt::xtensor<T, 3> _vals({X.shape(0), 1, value_size});
539
540 // Tabulate 1st order derivatives of shape functions at interpolation coords
541 xt::xtensor<double, 3> dphi = xt::view(
542 cmap.tabulate(1, X), xt::range(1, tdim + 1), xt::all(), xt::all(), 0);
543
544 const std::function<void(const xtl::span<T>&,
545 const xtl::span<const std::uint32_t>&,
546 std::int32_t, int)>
547 apply_inverse_transpose_dof_transformation
548 = element->get_dof_transformation_function<T>(true, true);
549
550 // Get interpolation operator
551 const xt::xtensor<double, 2>& Pi = element->interpolation_operator();
552
553 using U_t = xt::xview<decltype(reference_data)&, std::size_t,
554 xt::xall<std::size_t>, xt::xall<std::size_t>>;
555 using J_t = xt::xview<decltype(J)&, std::size_t, xt::xall<std::size_t>,
556 xt::xall<std::size_t>>;
557 auto pull_back_fn = element->map_fn<U_t, U_t, J_t, J_t>();
558 for (std::size_t c = 0; c < cells.size(); ++c)
559 {
560 const std::int32_t cell = cells[c];
561 auto x_dofs = x_dofmap.links(cell);
562 for (int i = 0; i < num_dofs_g; ++i)
563 {
564 const int pos = 3 * x_dofs[i];
565 for (int j = 0; j < gdim; ++j)
566 coordinate_dofs(i, j) = x_g[pos + j];
567 }
568
569 // Compute J, detJ and K
570 J.fill(0);
571 for (std::size_t p = 0; p < X.shape(0); ++p)
572 {
573 cmap.compute_jacobian(xt::view(dphi, xt::all(), p, xt::all()),
574 coordinate_dofs,
575 xt::view(J, p, xt::all(), xt::all()));
576 cmap.compute_jacobian_inverse(xt::view(J, p, xt::all(), xt::all()),
577 xt::view(K, p, xt::all(), xt::all()));
578 detJ[p] = cmap.compute_jacobian_determinant(
579 xt::view(J, p, xt::all(), xt::all()));
580 }
581
582 xtl::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
583 for (int k = 0; k < element_bs; ++k)
584 {
585 // Extract computed expression values for element block k
586 for (int m = 0; m < value_size; ++m)
587 {
588 std::copy_n(std::next(_f.begin(), f_shape1 * (k * value_size + m)
589 + c * X.shape(0)),
590 X.shape(0), xt::view(_vals, xt::all(), 0, m).begin());
591 }
592
593 // Get element degrees of freedom for block
594 for (std::size_t i = 0; i < X.shape(0); ++i)
595 {
596 auto _K = xt::view(K, i, xt::all(), xt::all());
597 auto _J = xt::view(J, i, xt::all(), xt::all());
598 auto _u = xt::view(_vals, i, xt::all(), xt::all());
599 auto _U = xt::view(reference_data, i, xt::all(), xt::all());
600 pull_back_fn(_U, _u, _K, 1.0 / detJ[i], _J);
601 }
602
603 auto ref_data = xt::view(reference_data, xt::all(), 0, xt::all());
604 impl::interpolation_apply(Pi, ref_data, _coeffs, element_bs);
605 apply_inverse_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
606
607 // Copy interpolation dofs into coefficient vector
608 assert(_coeffs.size() == num_scalar_dofs);
609 for (int i = 0; i < num_scalar_dofs; ++i)
610 {
611 const int dof = i * element_bs + k;
612 std::div_t pos = std::div(dof, dofmap_bs);
613 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
614 }
615 }
616 }
617 }
618}
619
625template <typename T>
627 const xtl::span<const std::int32_t>& cells)
628{
629 assert(u.function_space());
630 assert(v.function_space());
631 std::shared_ptr<const mesh::Mesh> mesh = u.function_space()->mesh();
632 assert(mesh);
633
634 auto cell_map0 = mesh->topology().index_map(mesh->topology().dim());
635 assert(cell_map0);
636 std::size_t num_cells0 = cell_map0->size_local() + cell_map0->num_ghosts();
637 if (u.function_space() == v.function_space() and cells.size() == num_cells0)
638 {
639 // Same function spaces and on whole mesh
640 xtl::span<T> u1_array = u.x()->mutable_array();
641 xtl::span<const T> u0_array = v.x()->array();
642 std::copy(u0_array.begin(), u0_array.end(), u1_array.begin());
643 }
644 else
645 {
646 // Get mesh and check that functions share the same mesh
647 if (mesh != v.function_space()->mesh())
648 {
649 throw std::runtime_error(
650 "Interpolation on different meshes not supported (yet).");
651 }
652
653 // Get elements and check value shape
654 auto element0 = v.function_space()->element();
655 assert(element0);
656 auto element1 = u.function_space()->element();
657 assert(element1);
658 if (element0->value_shape() != element1->value_shape())
659 {
660 throw std::runtime_error(
661 "Interpolation: elements have different value dimensions");
662 }
663
664 if (*element1 == *element0)
665 {
666 // Same element, different dofmaps (or just a subset of cells)
667
668 const int tdim = mesh->topology().dim();
669 auto cell_map = mesh->topology().index_map(tdim);
670 assert(cell_map);
671
672 assert(element1->block_size() == element0->block_size());
673
674 // Get dofmaps
675 std::shared_ptr<const fem::DofMap> dofmap0 = v.function_space()->dofmap();
676 assert(dofmap0);
677 std::shared_ptr<const fem::DofMap> dofmap1 = u.function_space()->dofmap();
678 assert(dofmap1);
679
680 xtl::span<T> u1_array = u.x()->mutable_array();
681 xtl::span<const T> u0_array = v.x()->array();
682
683 // Iterate over mesh and interpolate on each cell
684 const int bs0 = dofmap0->bs();
685 const int bs1 = dofmap1->bs();
686 for (auto c : cells)
687 {
688 xtl::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
689 xtl::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
690 assert(bs0 * dofs0.size() == bs1 * dofs1.size());
691 for (std::size_t i = 0; i < dofs0.size(); ++i)
692 {
693 for (int k = 0; k < bs0; ++k)
694 {
695 int index = bs0 * i + k;
696 std::div_t dv1 = std::div(index, bs1);
697 u1_array[bs1 * dofs1[dv1.quot] + dv1.rem]
698 = u0_array[bs0 * dofs0[i] + k];
699 }
700 }
701 }
702 }
703 else if (element1->map_type() == element0->map_type())
704 {
705 // Different elements, same basis function map type
706 impl::interpolate_same_map(u, v, cells);
707 }
708 else
709 {
710 // Different elements with different maps for basis functions
711 impl::interpolate_nonmatching_maps(u, v, cells);
712 }
713 }
714}
715
716} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
Definition: CoordinateElement.h:31
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:107
static double compute_jacobian_determinant(const U &J)
Compute the determinant of the Jacobian.
Definition: CoordinateElement.h:130
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition: CoordinateElement.h:116
int dim() const
The dimension of the geometry element space.
Definition: CoordinateElement.cpp:207
xt::xtensor< double, 4 > tabulate(int nd, const xt::xtensor< double, 2 > &X) const
Evaluate basis values and derivatives at set of points.
Definition: CoordinateElement.cpp:50
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:46
xtl::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:131
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
xt::xtensor< double, 2 > interpolation_coords(const fem::FiniteElement &element, const mesh::Mesh &mesh, const xtl::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:17
void interpolate(Function< T > &u, const xt::xarray< T > &f, const xtl::span< const std::int32_t > &cells)
Interpolate an expression f(x) in a finite element space.
Definition: interpolate.h:383