DOLFIN-X
DOLFIN-X C++ interface
interpolate.h
1 // Copyright (C) 2020-2021 Garth N. Wells
2 //
3 // This file is part of DOLFINX (https://www.fenicsproject.org)
4 //
5 // SPDX-License-Identifier: LGPL-3.0-or-later
6 
7 #pragma once
8 
9 #include "FunctionSpace.h"
10 #include <dolfinx/common/span.hpp>
11 #include <dolfinx/fem/DofMap.h>
12 #include <dolfinx/fem/FiniteElement.h>
13 #include <dolfinx/mesh/Mesh.h>
14 #include <functional>
15 #include <variant>
16 #include <xtensor/xtensor.hpp>
17 
18 namespace dolfinx::fem
19 {
20 
21 template <typename T>
22 class Function;
23 
34 xt::xtensor<double, 2>
35 interpolation_coords(const fem::FiniteElement& element, const mesh::Mesh& mesh,
36  const tcb::span<const std::int32_t>& cells);
37 
42 template <typename T>
43 void interpolate(Function<T>& u, const Function<T>& v);
44 
56 template <typename T>
57 void interpolate(
58  Function<T>& u,
59  const std::function<xt::xarray<T>(const xt::xtensor<double, 2>&)>& f,
60  const xt::xtensor<double, 2>& x,
61  const tcb::span<const std::int32_t>& cells);
62 
79 template <typename T>
80 void interpolate_c(
81  Function<T>& u,
82  const std::function<void(xt::xarray<T>&, const xt::xtensor<double, 2>&)>& f,
83  const xt::xtensor<double, 2>& x,
84  const tcb::span<const std::int32_t>& cells);
85 
86 namespace detail
87 {
88 
89 template <typename T>
90 void interpolate_from_any(Function<T>& u, const Function<T>& v)
91 {
92  assert(v.function_space());
93  const auto element = u.function_space()->element();
94  assert(element);
95  if (v.function_space()->element()->hash() != element->hash())
96  {
97  throw std::runtime_error("Restricting finite elements function in "
98  "different elements not supported.");
99  }
100 
101  const auto mesh = u.function_space()->mesh();
102  assert(mesh);
103  assert(v.function_space()->mesh());
104  if (mesh->id() != v.function_space()->mesh()->id())
105  {
106  throw std::runtime_error(
107  "Interpolation on different meshes not supported (yet).");
108  }
109  const int tdim = mesh->topology().dim();
110 
111  // Get dofmaps
112  assert(v.function_space());
113  std::shared_ptr<const fem::DofMap> dofmap_v = v.function_space()->dofmap();
114  assert(dofmap_v);
115  auto map = mesh->topology().index_map(tdim);
116  assert(map);
117 
118  std::vector<T>& coeffs = u.x()->mutable_array();
119 
120  // Iterate over mesh and interpolate on each cell
121  const auto dofmap_u = u.function_space()->dofmap();
122  const std::vector<T>& v_array = v.x()->array();
123  const int num_cells = map->size_local() + map->num_ghosts();
124  const int bs = dofmap_v->bs();
125  assert(bs == dofmap_u->bs());
126  for (int c = 0; c < num_cells; ++c)
127  {
128  tcb::span<const std::int32_t> dofs_v = dofmap_v->cell_dofs(c);
129  tcb::span<const std::int32_t> cell_dofs = dofmap_u->cell_dofs(c);
130  assert(dofs_v.size() == cell_dofs.size());
131  for (std::size_t i = 0; i < dofs_v.size(); ++i)
132  {
133  for (int k = 0; k < bs; ++k)
134  coeffs[bs * cell_dofs[i] + k] = v_array[bs * dofs_v[i] + k];
135  }
136  }
137 }
138 
139 } // namespace detail
140 
141 //----------------------------------------------------------------------------
142 template <typename T>
144 {
145  assert(u.function_space());
146  const std::shared_ptr<const FiniteElement> element
147  = u.function_space()->element();
148  assert(element);
149 
150  // Check that function ranks match
151  if (int rank_v = v.function_space()->element()->value_rank();
152  element->value_rank() != rank_v)
153  {
154  throw std::runtime_error("Cannot interpolate function into function space. "
155  "Rank of function ("
156  + std::to_string(rank_v)
157  + ") does not match rank of function space ("
158  + std::to_string(element->value_rank()) + ")");
159  }
160 
161  // Check that function dimension match
162  for (int i = 0; i < element->value_rank(); ++i)
163  {
164  if (int v_dim = v.function_space()->element()->value_dimension(i);
165  element->value_dimension(i) != v_dim)
166  {
167  throw std::runtime_error(
168  "Cannot interpolate function into function space. "
169  "Dimension "
170  + std::to_string(i) + " of function (" + std::to_string(v_dim)
171  + ") does not match dimension " + std::to_string(i)
172  + " of function space(" + std::to_string(element->value_dimension(i))
173  + ")");
174  }
175  }
176 
177  detail::interpolate_from_any(u, v);
178 }
179 //----------------------------------------------------------------------------
180 template <typename T>
182  Function<T>& u,
183  const std::function<xt::xarray<T>(const xt::xtensor<double, 2>&)>& f,
184  const xt::xtensor<double, 2>& x, const tcb::span<const std::int32_t>& cells)
185 {
186  const std::shared_ptr<const FiniteElement> element
187  = u.function_space()->element();
188  assert(element);
189  const int element_bs = element->block_size();
190  if (int num_sub = element->num_sub_elements();
191  num_sub > 0 and num_sub != element_bs)
192  {
193  throw std::runtime_error("Cannot directly interpolate a mixed space. "
194  "Interpolate into subspaces.");
195  }
196 
197  // Get mesh
198  assert(u.function_space());
199  auto mesh = u.function_space()->mesh();
200  assert(mesh);
201 
202  const int gdim = mesh->geometry().dim();
203  const int tdim = mesh->topology().dim();
204 
205  // Get the interpolation points on the reference cells
206  const array2d<double> X = element->interpolation_points();
207 
208  if (X.shape[0] == 0)
209  throw std::runtime_error(
210  "Interpolation into this space is not yet supported.");
211 
212  mesh->topology_mutable().create_entity_permutations();
213  const std::vector<std::uint32_t>& cell_info
214  = mesh->topology().get_cell_permutation_info();
215 
216  // Evaluate function at physical points. The returned array has a
217  // number of rows equal to the number of components of the function,
218  // and the number of columns is equal to the number of evaluation
219  // points.
220 
221  xt::xarray<T> values = f(x);
222 
223  if (values.dimension() == 1)
224  {
225  if (element->value_size() != 1)
226  throw std::runtime_error("Interpolation data has the wrong shape.");
227  values.reshape({element->value_size(), x.shape(1)});
228  }
229 
230  if (values.shape(0) != element->value_size())
231  throw std::runtime_error("Interpolation data has the wrong shape.");
232 
233  if (values.shape(1) != cells.size() * X.shape[0])
234  throw std::runtime_error("Interpolation data has the wrong shape.");
235 
236  // Get dofmap
237  const auto dofmap = u.function_space()->dofmap();
238  assert(dofmap);
239  const int dofmap_bs = dofmap->bs();
240 
241  // Loop over cells and compute interpolation dofs
242  const int num_scalar_dofs = element->space_dimension() / element_bs;
243  const int value_size = element->value_size() / element_bs;
244 
245  std::vector<T>& coeffs = u.x()->mutable_array();
246  std::vector<T> _coeffs(num_scalar_dofs);
247 
248  // This assumes that any element with an identity interpolation matrix is a
249  // point evaluation
250  if (element->interpolation_ident())
251  {
252  for (std::int32_t c : cells)
253  {
254  tcb::span<const std::int32_t> dofs = dofmap->cell_dofs(c);
255  for (int k = 0; k < element_bs; ++k)
256  {
257  for (int i = 0; i < num_scalar_dofs; ++i)
258  _coeffs[i] = values(k, c * num_scalar_dofs + i);
259  element->apply_inverse_transpose_dof_transformation(_coeffs.data(),
260  cell_info[c], 1);
261  for (int i = 0; i < num_scalar_dofs; ++i)
262  {
263  const int dof = i * element_bs + k;
264  std::div_t pos = std::div(dof, dofmap_bs);
265  coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
266  }
267  }
268  }
269  }
270  else
271  {
272  // Get coordinate map
273  const fem::CoordinateElement& cmap = mesh->geometry().cmap();
274 
275  // Get geometry data
276  const graph::AdjacencyList<std::int32_t>& x_dofmap
277  = mesh->geometry().dofmap();
278  // FIXME: Add proper interface for num coordinate dofs
279  const int num_dofs_g = x_dofmap.num_links(0);
280  const array2d<double>& x_g = mesh->geometry().x();
281 
282  // Create data structures for Jacobian info
283  array2d<double> x_cell(X.shape[0], gdim);
284  std::vector<double> J(X.shape[0] * gdim * tdim);
285  std::vector<double> detJ(X.shape[0]);
286  std::vector<double> K(X.shape[0] * tdim * gdim);
287  array2d<double> X_ref(X.shape[0], tdim);
288 
289  array2d<double> coordinate_dofs(num_dofs_g, gdim);
290 
291  array2d<T> reference_data(value_size, X.shape[0]);
292  array2d<T> _vals(value_size, X.shape[0]);
293 
294  // Tabulate 0th and 1st order derivatives of shape functions at
295  // interpolation coords
296  xt::xtensor<double, 4> tabulated_data = cmap.tabulate(1, X);
297  for (std::int32_t c : cells)
298  {
299  auto x_dofs = x_dofmap.links(c);
300  for (int i = 0; i < num_dofs_g; ++i)
301  for (int j = 0; j < gdim; ++j)
302  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
303 
304  // Compute J, detJ and K
305  cmap.compute_jacobian_data(tabulated_data, X, coordinate_dofs, J, detJ,
306  K);
307 
308  tcb::span<const std::int32_t> dofs = dofmap->cell_dofs(c);
309  for (int k = 0; k < element_bs; ++k)
310  {
311  // Extract computed expression values for element block k
312  for (int m = 0; m < value_size; ++m)
313  {
314  std::copy_n(&values(k * value_size + m, c * X.shape[0]), X.shape[0],
315  _vals.row(m).begin());
316  }
317 
318  // Get element degrees of freedom for block
319  element->map_pull_back(reference_data.data(), _vals.data(), J.data(),
320  detJ.data(), K.data(), gdim, value_size, 1,
321  X.shape[0]);
322 
323  element->interpolate(reference_data, tcb::make_span(_coeffs));
324  element->apply_inverse_transpose_dof_transformation(_coeffs.data(),
325  cell_info[c], 1);
326 
327  assert(_coeffs.size() == num_scalar_dofs);
328 
329  // Copy interpolation dofs into coefficient vector
330  for (int i = 0; i < num_scalar_dofs; ++i)
331  {
332  const int dof = i * element_bs + k;
333  std::div_t pos = std::div(dof, dofmap_bs);
334  coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
335  }
336  }
337  }
338  }
339 }
340 //----------------------------------------------------------------------------
341 template <typename T>
343  Function<T>& u,
344  const std::function<void(xt::xarray<T>&, const xt::xtensor<double, 2>&)>& f,
345  const xt::xtensor<double, 2>& x, const tcb::span<const std::int32_t>& cells)
346 {
347  const std::shared_ptr<const FiniteElement> element
348  = u.function_space()->element();
349  assert(element);
350  std::vector<int> vshape(element->value_rank(), 1);
351  for (std::size_t i = 0; i < vshape.size(); ++i)
352  vshape[i] = element->value_dimension(i);
353  const std::size_t value_size = std::accumulate(
354  std::begin(vshape), std::end(vshape), 1, std::multiplies<>());
355 
356  auto fn = [value_size, &f](const xt::xtensor<double, 2>& x) {
357  xt::xarray<T> values = xt::empty<T>({value_size, x.shape(1)});
358  f(values, x);
359  return values;
360  };
361 
362  interpolate<T>(u, fn, x, cells);
363 }
364 //----------------------------------------------------------------------------
365 
366 } // namespace dolfinx::fem
This class provides a dynamic 2-dimensional row-wise array data structure.
Definition: array2d.h:21
std::array< size_type, 2 > shape
The shape of the array.
Definition: array2d.h:157
constexpr tcb::span< value_type > row(size_type i)
Access a row in the array.
Definition: array2d.h:116
constexpr value_type * data() noexcept
Get pointer to the first element of the underlying storage.
Definition: array2d.h:133
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:31
xt::xtensor< double, 4 > tabulate(int n, const array2d< double > &X) const
Tabulate shape functions up to n-th order derivative at points X in the reference geometry Note: Dyna...
Definition: CoordinateElement.cpp:248
void compute_jacobian_data(const xt::xtensor< double, 4 > &tabulated_data, const array2d< double > &X, const array2d< double > &cell_geometry, std::vector< double > &J, tcb::span< double > detJ, std::vector< double > &K) const
Compute J, K and detJ for a cell with given geometry, and the basis functions and first order derivat...
Definition: CoordinateElement.cpp:255
This class represents a function in a finite element function space , given by.
Definition: Function.h:46
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:154
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:195
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:68
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:151
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:141
Finite element method functionality.
Definition: assemble_matrix_impl.h:23
xt::xtensor< double, 2 > interpolation_coords(const fem::FiniteElement &element, const mesh::Mesh &mesh, const tcb::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_c(Function< T > &u, const std::function< void(xt::xarray< T > &, const xt::xtensor< double, 2 > &)> &f, const xt::xtensor< double, 2 > &x, const tcb::span< const std::int32_t > &cells)
Interpolate an expression f(x)
Definition: interpolate.h:342
void interpolate(Function< T > &u, const Function< T > &v)
Interpolate a finite element Function (on possibly non-matching meshes) in another finite element spa...
Definition: interpolate.h:143