Note: this is documentation for an old release. View the latest documentation at docs.fenicsproject.org/dolfinx/v0.9.0/cpp/doxygen/df/dd3/interpolate_8h_source.html
DOLFINx  0.5.0
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 
21 namespace dolfinx::fem
22 {
23 template <typename T>
24 class Function;
25 
26 namespace impl
27 {
36 template <typename U, typename V, typename T>
37 void 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 
83 template <typename T>
84 void 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 
168 template <typename T>
169 void 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 
404 std::vector<double> interpolation_coords(const FiniteElement& element,
405  const mesh::Mesh& mesh,
406  std::span<const std::int32_t> cells);
407 
420 template <typename T>
421 void 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->interpolation_ident())
480  {
481  // Point evaluation element *and* the geometric map is the identity,
482  // e.g. not Piola mapped
483 
484  if (!element->map_ident())
485  throw std::runtime_error("Element does not have identity map.");
486 
487  auto apply_inv_transpose_dof_transformation
488  = element->get_dof_transformation_function<T>(true, true, true);
489 
490  // Loop over cells
491  for (std::size_t c = 0; c < cells.size(); ++c)
492  {
493  const std::int32_t cell = cells[c];
494  std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
495  for (int k = 0; k < element_bs; ++k)
496  {
497  // num_scalar_dofs is the number of interpolation points per
498  // cell in this case (interpolation matrix is identity)
499  std::copy_n(std::next(f.begin(), k * f_shape1 + c * num_scalar_dofs),
500  num_scalar_dofs, _coeffs.begin());
501  apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
502  for (int i = 0; i < num_scalar_dofs; ++i)
503  {
504  const int dof = i * element_bs + k;
505  std::div_t pos = std::div(dof, dofmap_bs);
506  coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
507  }
508  }
509  }
510  }
511  else if (element->map_ident())
512  {
513  // Not a point evaluation, but the geometric map is the identity,
514  // e.g. not Piola mapped
515 
516  if (_f.extent(0) != 1)
517  throw std::runtime_error("Interpolation data has the wrong shape.");
518 
519  // Get interpolation operator
520  const auto [_Pi, pi_shape] = element->interpolation_operator();
521  cmdspan2_t Pi(_Pi.data(), pi_shape);
522  const std::size_t num_interp_points = Pi.extent(1);
523  assert(Pi.extent(0) == num_scalar_dofs);
524 
525  auto apply_inv_transpose_dof_transformation
526  = element->get_dof_transformation_function<T>(true, true, true);
527 
528  // Loop over cells
529  std::vector<T> ref_data_b(num_interp_points);
530  stdex::mdspan<T, stdex::extents<std::size_t, stdex::dynamic_extent, 1>>
531  ref_data(ref_data_b.data(), num_interp_points, 1);
532  for (std::size_t c = 0; c < cells.size(); ++c)
533  {
534  const std::int32_t cell = cells[c];
535  std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
536  for (int k = 0; k < element_bs; ++k)
537  {
538  std::copy_n(std::next(f.begin(), k * f_shape1 + c * num_interp_points),
539  num_interp_points, ref_data_b.begin());
540  impl::interpolation_apply(Pi, ref_data, std::span(_coeffs), element_bs);
541  apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
542  for (int i = 0; i < num_scalar_dofs; ++i)
543  {
544  const int dof = i * element_bs + k;
545  std::div_t pos = std::div(dof, dofmap_bs);
546  coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
547  }
548  }
549  }
550  }
551  else
552  {
553  // Get the interpolation points on the reference cells
554  const auto [X, Xshape] = element->interpolation_points();
555  if (X.empty())
556  {
557  throw std::runtime_error(
558  "Interpolation into this space is not yet supported.");
559  }
560 
561  if (_f.extent(1) != cells.size() * Xshape[0])
562  throw std::runtime_error("Interpolation data has the wrong shape.");
563 
564  // Get coordinate map
565  const CoordinateElement& cmap = mesh->geometry().cmap();
566 
567  // Get geometry data
568  const graph::AdjacencyList<std::int32_t>& x_dofmap
569  = mesh->geometry().dofmap();
570  const int num_dofs_g = cmap.dim();
571  std::span<const double> x_g = mesh->geometry().x();
572 
573  // Create data structures for Jacobian info
574  std::vector<double> J_b(Xshape[0] * gdim * tdim);
575  mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
576  std::vector<double> K_b(Xshape[0] * tdim * gdim);
577  mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
578  std::vector<double> detJ(Xshape[0]);
579  std::vector<double> det_scratch(2 * gdim * tdim);
580 
581  std::vector<double> coord_dofs_b(num_dofs_g * gdim);
582  mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
583 
584  std::vector<T> ref_data_b(Xshape[0] * 1 * value_size);
585  stdex::mdspan<T, stdex::dextents<std::size_t, 3>> ref_data(
586  ref_data_b.data(), Xshape[0], 1, value_size);
587 
588  std::vector<T> _vals_b(Xshape[0] * 1 * value_size);
589  stdex::mdspan<T, stdex::dextents<std::size_t, 3>> _vals(
590  _vals_b.data(), Xshape[0], 1, value_size);
591 
592  // Tabulate 1st derivative of shape functions at interpolation
593  // coords
594  std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(1, Xshape[0]);
595  std::vector<double> phi_b(
596  std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
597  cmdspan4_t phi(phi_b.data(), phi_shape);
598  cmap.tabulate(1, X, Xshape, phi_b);
599  auto dphi = stdex::submdspan(phi, std::pair(1, tdim + 1),
600  stdex::full_extent, stdex::full_extent, 0);
601 
602  const std::function<void(const std::span<T>&,
603  const std::span<const std::uint32_t>&,
604  std::int32_t, int)>
605  apply_inverse_transpose_dof_transformation
606  = element->get_dof_transformation_function<T>(true, true);
607 
608  // Get interpolation operator
609  const auto [_Pi, pi_shape] = element->interpolation_operator();
610  cmdspan2_t Pi(_Pi.data(), pi_shape);
611 
612  namespace stdex = std::experimental;
613  using u_t = stdex::mdspan<const T, stdex::dextents<std::size_t, 2>>;
614  using U_t = stdex::mdspan<T, stdex::dextents<std::size_t, 2>>;
615  using J_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
616  using K_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
617  auto pull_back_fn = element->basix_element().map_fn<U_t, u_t, J_t, K_t>();
618 
619  for (std::size_t c = 0; c < cells.size(); ++c)
620  {
621  const std::int32_t cell = cells[c];
622  auto x_dofs = x_dofmap.links(cell);
623  for (int i = 0; i < num_dofs_g; ++i)
624  {
625  const int pos = 3 * x_dofs[i];
626  for (int j = 0; j < gdim; ++j)
627  coord_dofs(i, j) = x_g[pos + j];
628  }
629 
630  // Compute J, detJ and K
631  std::fill(J_b.begin(), J_b.end(), 0);
632  for (std::size_t p = 0; p < Xshape[0]; ++p)
633  {
634  auto _dphi
635  = stdex::submdspan(dphi, stdex::full_extent, p, stdex::full_extent);
636  auto _J
637  = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
638  cmap.compute_jacobian(_dphi, coord_dofs, _J);
639  auto _K
640  = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
641  cmap.compute_jacobian_inverse(_J, _K);
642  detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
643  }
644 
645  std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
646  for (int k = 0; k < element_bs; ++k)
647  {
648  // Extract computed expression values for element block k
649  for (int m = 0; m < value_size; ++m)
650  {
651  for (std::size_t k0 = 0; k0 < Xshape[0]; ++k0)
652  {
653  _vals(k0, 0, m)
654  = f[f_shape1 * (k * value_size + m) + c * Xshape[0] + k0];
655  }
656  }
657 
658  // Get element degrees of freedom for block
659  for (std::size_t i = 0; i < Xshape[0]; ++i)
660  {
661  auto _u = stdex::submdspan(_vals, i, stdex::full_extent,
662  stdex::full_extent);
663  auto _U = stdex::submdspan(ref_data, i, stdex::full_extent,
664  stdex::full_extent);
665  auto _K
666  = stdex::submdspan(K, i, stdex::full_extent, stdex::full_extent);
667  auto _J
668  = stdex::submdspan(J, i, stdex::full_extent, stdex::full_extent);
669  pull_back_fn(_U, _u, _K, 1.0 / detJ[i], _J);
670  }
671 
672  auto ref = stdex::submdspan(ref_data, stdex::full_extent, 0,
673  stdex::full_extent);
674  impl::interpolation_apply(Pi, ref, std::span(_coeffs), element_bs);
675  apply_inverse_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
676 
677  // Copy interpolation dofs into coefficient vector
678  assert(_coeffs.size() == num_scalar_dofs);
679  for (int i = 0; i < num_scalar_dofs; ++i)
680  {
681  const int dof = i * element_bs + k;
682  std::div_t pos = std::div(dof, dofmap_bs);
683  coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
684  }
685  }
686  }
687  }
688 }
689 
695 template <typename T>
697  const std::span<const std::int32_t>& cells)
698 {
699  assert(u.function_space());
700  assert(v.function_space());
701  std::shared_ptr<const mesh::Mesh> mesh = u.function_space()->mesh();
702  assert(mesh);
703 
704  auto cell_map0 = mesh->topology().index_map(mesh->topology().dim());
705  assert(cell_map0);
706  std::size_t num_cells0 = cell_map0->size_local() + cell_map0->num_ghosts();
707  if (u.function_space() == v.function_space() and cells.size() == num_cells0)
708  {
709  // Same function spaces and on whole mesh
710  std::span<T> u1_array = u.x()->mutable_array();
711  std::span<const T> u0_array = v.x()->array();
712  std::copy(u0_array.begin(), u0_array.end(), u1_array.begin());
713  }
714  else
715  {
716  // Get mesh and check that functions share the same mesh
717  if (mesh != v.function_space()->mesh())
718  {
719  throw std::runtime_error(
720  "Interpolation on different meshes not supported (yet).");
721  }
722 
723  // Get elements and check value shape
724  auto element0 = v.function_space()->element();
725  assert(element0);
726  auto element1 = u.function_space()->element();
727  assert(element1);
728  if (element0->value_shape().size() != element1->value_shape().size()
729  or !std::equal(element0->value_shape().begin(),
730  element0->value_shape().end(),
731  element1->value_shape().begin()))
732  {
733  throw std::runtime_error(
734  "Interpolation: elements have different value dimensions");
735  }
736 
737  if (*element1 == *element0)
738  {
739  // Same element, different dofmaps (or just a subset of cells)
740 
741  const int tdim = mesh->topology().dim();
742  auto cell_map = mesh->topology().index_map(tdim);
743  assert(cell_map);
744 
745  assert(element1->block_size() == element0->block_size());
746 
747  // Get dofmaps
748  std::shared_ptr<const DofMap> dofmap0 = v.function_space()->dofmap();
749  assert(dofmap0);
750  std::shared_ptr<const DofMap> dofmap1 = u.function_space()->dofmap();
751  assert(dofmap1);
752 
753  std::span<T> u1_array = u.x()->mutable_array();
754  std::span<const T> u0_array = v.x()->array();
755 
756  // Iterate over mesh and interpolate on each cell
757  const int bs0 = dofmap0->bs();
758  const int bs1 = dofmap1->bs();
759  for (auto c : cells)
760  {
761  std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
762  std::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(c);
763  assert(bs0 * dofs0.size() == bs1 * dofs1.size());
764  for (std::size_t i = 0; i < dofs0.size(); ++i)
765  {
766  for (int k = 0; k < bs0; ++k)
767  {
768  int index = bs0 * i + k;
769  std::div_t dv1 = std::div(index, bs1);
770  u1_array[bs1 * dofs1[dv1.quot] + dv1.rem]
771  = u0_array[bs0 * dofs0[i] + k];
772  }
773  }
774  }
775  }
776  else if (element1->map_type() == element0->map_type())
777  {
778  // Different elements, same basis function map type
779  impl::interpolate_same_map(u, v, cells);
780  }
781  else
782  {
783  // Different elements with different maps for basis functions
784  impl::interpolate_nonmatching_maps(u, v, cells);
785  }
786  }
787 }
788 
789 } // 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 la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:142
std::shared_ptr< const FunctionSpace > function_space() const
Access the function space.
Definition: Function.h:136
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