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 <Eigen/Core>
11 #include <dolfinx/common/span.hpp>
12 #include <dolfinx/fem/DofMap.h>
13 #include <dolfinx/fem/FiniteElement.h>
14 #include <dolfinx/mesh/Mesh.h>
15 #include <functional>
16 #include <variant>
17 
18 namespace dolfinx::fem
19 {
20 
21 template <typename T>
22 class Function;
23 
34 array2d<double>
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(Function<T>& u,
58  const std::function<std::variant<std::vector<T>, array2d<T>>(
59  const array2d<double>&)>& f,
60  const array2d<double>& 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(array2d<T>&, const array2d<double>&)>& f,
83  const array2d<double>& x, const tcb::span<const std::int32_t>& cells);
84 
85 namespace detail
86 {
87 
88 template <typename T>
89 void interpolate_from_any(Function<T>& u, const Function<T>& v)
90 {
91  assert(v.function_space());
92  const auto element = u.function_space()->element();
93  assert(element);
94  if (v.function_space()->element()->hash() != element->hash())
95  {
96  throw std::runtime_error("Restricting finite elements function in "
97  "different elements not supported.");
98  }
99 
100  const auto mesh = u.function_space()->mesh();
101  assert(mesh);
102  assert(v.function_space()->mesh());
103  if (mesh->id() != v.function_space()->mesh()->id())
104  {
105  throw std::runtime_error(
106  "Interpolation on different meshes not supported (yet).");
107  }
108  const int tdim = mesh->topology().dim();
109 
110  // Get dofmaps
111  assert(v.function_space());
112  std::shared_ptr<const fem::DofMap> dofmap_v = v.function_space()->dofmap();
113  assert(dofmap_v);
114  auto map = mesh->topology().index_map(tdim);
115  assert(map);
116 
117  std::vector<T>& coeffs = u.x()->mutable_array();
118 
119  // Iterate over mesh and interpolate on each cell
120  const auto dofmap_u = u.function_space()->dofmap();
121  const std::vector<T>& v_array = v.x()->array();
122  const int num_cells = map->size_local() + map->num_ghosts();
123  const int bs = dofmap_v->bs();
124  assert(bs == dofmap_u->bs());
125  for (int c = 0; c < num_cells; ++c)
126  {
127  tcb::span<const std::int32_t> dofs_v = dofmap_v->cell_dofs(c);
128  tcb::span<const std::int32_t> cell_dofs = dofmap_u->cell_dofs(c);
129  assert(dofs_v.size() == cell_dofs.size());
130  for (std::size_t i = 0; i < dofs_v.size(); ++i)
131  {
132  for (int k = 0; k < bs; ++k)
133  coeffs[bs * cell_dofs[i] + k] = v_array[bs * dofs_v[i] + k];
134  }
135  }
136 }
137 
138 } // namespace detail
139 
140 //----------------------------------------------------------------------------
141 template <typename T>
143 {
144  assert(u.function_space());
145  const auto element = u.function_space()->element();
146  assert(element);
147 
148  // Check that function ranks match
149  if (int rank_v = v.function_space()->element()->value_rank();
150  element->value_rank() != rank_v)
151  {
152  throw std::runtime_error("Cannot interpolate function into function space. "
153  "Rank of function ("
154  + std::to_string(rank_v)
155  + ") does not match rank of function space ("
156  + std::to_string(element->value_rank()) + ")");
157  }
158 
159  // Check that function dimension match
160  for (int i = 0; i < element->value_rank(); ++i)
161  {
162  if (int v_dim = v.function_space()->element()->value_dimension(i);
163  element->value_dimension(i) != v_dim)
164  {
165  throw std::runtime_error(
166  "Cannot interpolate function into function space. "
167  "Dimension "
168  + std::to_string(i) + " of function (" + std::to_string(v_dim)
169  + ") does not match dimension " + std::to_string(i)
170  + " of function space(" + std::to_string(element->value_dimension(i))
171  + ")");
172  }
173  }
174 
175  detail::interpolate_from_any(u, v);
176 }
177 //----------------------------------------------------------------------------
178 template <typename T>
180  const std::function<std::variant<std::vector<T>, array2d<T>>(
181  const array2d<double>&)>& f,
182  const array2d<double>& x,
183  const tcb::span<const std::int32_t>& cells)
184 {
185  const auto element = u.function_space()->element();
186  assert(element);
187  const int element_bs = element->block_size();
188  if (int num_sub = element->num_sub_elements();
189  num_sub > 0 and num_sub != element_bs)
190  {
191  throw std::runtime_error("Cannot directly interpolate a mixed space. "
192  "Interpolate into subspaces.");
193  }
194 
195  // Get mesh
196  assert(u.function_space());
197  auto mesh = u.function_space()->mesh();
198  assert(mesh);
199 
200  const int gdim = mesh->geometry().dim();
201  const int tdim = mesh->topology().dim();
202 
203  // Get the interpolation points on the reference cells
204  const array2d<double> X = element->interpolation_points();
205 
206  if (X.shape[0] == 0)
207  throw std::runtime_error(
208  "Interpolation into this space is not yet supported.");
209 
210  mesh->topology_mutable().create_entity_permutations();
211  const std::vector<std::uint32_t>& cell_info
212  = mesh->topology().get_cell_permutation_info();
213 
214  // Evaluate function at physical points. The returned array has a
215  // number of rows equal to the number of components of the function,
216  // and the number of columns is equal to the number of evaluation
217  // points.
218 
219  // TODO: Copies and memory allocation could be avoided with a 'span2d'
220  // class, or by just pointing to the data
221  array2d<T> values(element->value_size(), x.shape[1]);
222  std::variant<std::vector<T>, array2d<T>> values_v = f(x);
223  if (std::holds_alternative<array2d<T>>(values_v))
224  {
225  values = std::get<1>(values_v);
226  if (values.shape[0] != element->value_size())
227  throw std::runtime_error("Interpolation data has the wrong shape.");
228  }
229  else
230  {
231  if (element->value_size() != 1)
232  throw std::runtime_error("Interpolation data has the wrong shape.");
233  std::copy(std::get<0>(values_v).begin(), std::get<0>(values_v).end(),
234  values.data());
235  }
236 
237  if (values.shape[1] != cells.size() * X.shape[0])
238  throw std::runtime_error("Interpolation data has the wrong shape.");
239 
240  // Get dofmap
241  const auto dofmap = u.function_space()->dofmap();
242  assert(dofmap);
243  const int dofmap_bs = dofmap->bs();
244 
245  // Get coordinate map
246  const fem::CoordinateElement& cmap = mesh->geometry().cmap();
247 
248  // Get geometry data
249  const graph::AdjacencyList<std::int32_t>& x_dofmap
250  = mesh->geometry().dofmap();
251  // FIXME: Add proper interface for num coordinate dofs
252  const int num_dofs_g = x_dofmap.num_links(0);
253  const array2d<double>& x_g = mesh->geometry().x();
254 
255  // NOTE: The below loop over cells could be skipped for some elements,
256  // e.g. Lagrange, where the interpolation is just the identity
257 
258  // Loop over cells and compute interpolation dofs
259  const int num_scalar_dofs = element->space_dimension() / element_bs;
260  const int value_size = element->value_size() / element_bs;
261 
262  array2d<double> x_cell(X.shape[0], gdim);
263  std::vector<double> J(X.shape[0] * gdim * tdim);
264  std::vector<double> detJ(X.shape[0]);
265  std::vector<double> K(X.shape[0] * tdim * gdim);
266  array2d<double> X_ref(X.shape[0], tdim);
267 
268  array2d<double> coordinate_dofs(num_dofs_g, gdim);
269 
270  array2d<T> reference_data(value_size, X.shape[0]);
271 
272  std::vector<T>& coeffs = u.x()->mutable_array();
273  std::vector<T> _coeffs(num_scalar_dofs);
274  array2d<T> _vals(value_size, X.shape[0]);
275  for (std::int32_t c : cells)
276  {
277  auto x_dofs = x_dofmap.links(c);
278  for (int i = 0; i < num_dofs_g; ++i)
279  for (int j = 0; j < gdim; ++j)
280  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
281  cmap.push_forward(x_cell, X, coordinate_dofs);
282 
283  cmap.compute_reference_geometry(X_ref, J, detJ, K, x_cell, coordinate_dofs);
284 
285  auto dofs = dofmap->cell_dofs(c);
286 
287  for (int k = 0; k < element_bs; ++k)
288  {
289  // Extract computed expression values for element block k
290  for (int m = 0; m < value_size; ++m)
291  {
292  std::copy_n(&values(k * value_size + m, c * X.shape[0]), X.shape[0],
293  _vals.row(m).begin());
294  }
295 
296  // Get element degrees of freedom for block
297  element->map_pull_back(reference_data.data(), _vals.data(), J.data(),
298  detJ.data(), K.data(), gdim, value_size, 1,
299  X.shape[0]);
300 
301  element->interpolate(reference_data, tcb::make_span(_coeffs));
302  element->apply_inverse_transpose_dof_transformation(_coeffs.data(),
303  cell_info[c], 1);
304 
305  assert(_coeffs.size() == num_scalar_dofs);
306 
307  // Copy interpolation dofs into coefficient vector
308  for (int i = 0; i < num_scalar_dofs; ++i)
309  {
310  const int dof = i * element_bs + k;
311  std::div_t pos = std::div(dof, dofmap_bs);
312  coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
313  }
314  }
315  }
316 }
317 //----------------------------------------------------------------------------
318 template <typename T>
320  Function<T>& u,
321  const std::function<void(array2d<T>&, const array2d<double>&)>& f,
322  const array2d<double>& x, const tcb::span<const std::int32_t>& cells)
323 {
324  const auto element = u.function_space()->element();
325  assert(element);
326  std::vector<int> vshape(element->value_rank(), 1);
327  for (std::size_t i = 0; i < vshape.size(); ++i)
328  vshape[i] = element->value_dimension(i);
329  const int value_size = std::accumulate(std::begin(vshape), std::end(vshape),
330  1, std::multiplies<>());
331 
332  auto fn = [value_size, &f](const array2d<double>& x) {
333  array2d<T> values(value_size, x.shape[1]);
334  f(values, x);
335  return values;
336  };
337 
338  interpolate<T>(u, fn, x, cells);
339 }
340 //----------------------------------------------------------------------------
341 
342 } // 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:25
void push_forward(array2d< double > &x, const array2d< double > &X, const array2d< double > &cell_geometry) const
Compute physical coordinates x for points X in the reference configuration.
Definition: CoordinateElement.cpp:67
void compute_reference_geometry(array2d< double > &X, std::vector< double > &J, tcb::span< double > detJ, std::vector< double > &K, const array2d< double > &x, const array2d< double > &cell_geometry) const
Compute reference coordinates X, and J, detJ and K for physical coordinates x.
Definition: CoordinateElement.cpp:89
This class represents a function in a finite element function space , given by.
Definition: Function.h:44
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:152
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:193
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
array2d< double > 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:15
void interpolate_c(Function< T > &u, const std::function< void(array2d< T > &, const array2d< double > &)> &f, const array2d< double > &x, const tcb::span< const std::int32_t > &cells)
Interpolate an expression f(x)
Definition: interpolate.h:319
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:142