DOLFIN-X
DOLFIN-X C++ interface
interpolate.h
1 // Copyright (C) 2020 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 "Function.h"
10 #include "FunctionSpace.h"
11 #include <Eigen/Dense>
12 #include <dolfinx/fem/DofMap.h>
13 #include <dolfinx/fem/FiniteElement.h>
14 #include <dolfinx/mesh/Mesh.h>
15 #include <functional>
16 
17 namespace dolfinx::fem
18 {
19 
20 template <typename T>
21 class Function;
22 
27 template <typename T>
28 void interpolate(Function<T>& u, const Function<T>& v);
29 
33 template <typename T>
34 void interpolate(
35  Function<T>& u,
36  const std::function<
37  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
38  const Eigen::Ref<const Eigen::Array<double, 3, Eigen::Dynamic,
39  Eigen::RowMajor>>&)>& f);
40 
52 template <typename T>
53 void interpolate_c(
54  Function<T>& u,
55  const std::function<void(
56  Eigen::Ref<
57  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
58  const Eigen::Ref<const Eigen::Array<double, Eigen::Dynamic, 3,
59  Eigen::RowMajor>>&)>& f);
60 
61 namespace detail
62 {
63 
64 template <typename T>
65 void interpolate_from_any(Function<T>& u, const Function<T>& v)
66 {
67  assert(v.function_space());
68  const auto element = u.function_space()->element();
69  assert(element);
70  if (v.function_space()->element()->hash() != element->hash())
71  {
72  throw std::runtime_error("Restricting finite elements function in "
73  "different elements not supported.");
74  }
75 
76  const auto mesh = u.function_space()->mesh();
77  assert(mesh);
78  assert(v.function_space()->mesh());
79  if (mesh->id() != v.function_space()->mesh()->id())
80  {
81  throw std::runtime_error(
82  "Interpolation on different meshes not supported (yet).");
83  }
84  const int tdim = mesh->topology().dim();
85 
86  // Get dofmaps
87  assert(v.function_space());
88  std::shared_ptr<const fem::DofMap> dofmap_v = v.function_space()->dofmap();
89  assert(dofmap_v);
90  auto map = mesh->topology().index_map(tdim);
91  assert(map);
92 
93  std::vector<T>& coeffs = u.x()->mutable_array();
94 
95  // Iterate over mesh and interpolate on each cell
96  const auto dofmap_u = u.function_space()->dofmap();
97  const std::vector<T>& v_array = v.x()->array();
98  const int num_cells = map->size_local() + map->num_ghosts();
99  const int bs = dofmap_v->bs();
100  assert(bs == dofmap_u->bs());
101  for (int c = 0; c < num_cells; ++c)
102  {
103  tcb::span<const std::int32_t> dofs_v = dofmap_v->cell_dofs(c);
104  tcb::span<const std::int32_t> cell_dofs = dofmap_u->cell_dofs(c);
105  assert(dofs_v.size() == cell_dofs.size());
106  for (std::size_t i = 0; i < dofs_v.size(); ++i)
107  {
108  for (int k = 0; k < bs; ++k)
109  coeffs[bs * cell_dofs[i] + k] = v_array[bs * dofs_v[i] + k];
110  }
111  }
112 }
113 
114 } // namespace detail
115 
116 //----------------------------------------------------------------------------
117 template <typename T>
119 {
120  assert(u.function_space());
121  const auto element = u.function_space()->element();
122  assert(element);
123 
124  // Check that function ranks match
125  if (int rank_v = v.function_space()->element()->value_rank();
126  element->value_rank() != rank_v)
127  {
128  throw std::runtime_error("Cannot interpolate function into function space. "
129  "Rank of function ("
130  + std::to_string(rank_v)
131  + ") does not match rank of function space ("
132  + std::to_string(element->value_rank()) + ")");
133  }
134 
135  // Check that function dimension match
136  for (int i = 0; i < element->value_rank(); ++i)
137  {
138  if (int v_dim = v.function_space()->element()->value_dimension(i);
139  element->value_dimension(i) != v_dim)
140  {
141  throw std::runtime_error(
142  "Cannot interpolate function into function space. "
143  "Dimension "
144  + std::to_string(i) + " of function (" + std::to_string(v_dim)
145  + ") does not match dimension " + std::to_string(i)
146  + " of function space(" + std::to_string(element->value_dimension(i))
147  + ")");
148  }
149  }
150 
151  detail::interpolate_from_any(u, v);
152 }
153 //----------------------------------------------------------------------------
154 template <typename T>
156  Function<T>& u,
157  const std::function<
158  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
159  const Eigen::Ref<const Eigen::Array<double, 3, Eigen::Dynamic,
160  Eigen::RowMajor>>&)>& f)
161 {
162  using EigenMatrixRowXd
163  = Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
164 
165  const auto element = u.function_space()->element();
166  assert(element);
167  const int element_bs = element->block_size();
168 
169  if (int num_sub = element->num_sub_elements();
170  num_sub > 0 and num_sub != element_bs)
171  {
172  throw std::runtime_error("Cannot directly interpolate a mixed space. "
173  "Interpolate into subspaces.");
174  }
175 
176  // Get mesh
177  assert(u.function_space());
178  auto mesh = u.function_space()->mesh();
179  assert(mesh);
180  const int tdim = mesh->topology().dim();
181  const int gdim = mesh->geometry().dim();
182  auto cell_map = mesh->topology().index_map(tdim);
183  assert(cell_map);
184  const int num_cells = cell_map->size_local() + cell_map->num_ghosts();
185 
186  // Get mesh geometry data and the element coordinate map
187  const graph::AdjacencyList<std::int32_t>& x_dofmap
188  = mesh->geometry().dofmap();
189  const int num_dofs_g = x_dofmap.num_links(0);
190  const EigenMatrixRowXd& x_g = mesh->geometry().x();
191  const fem::CoordinateElement& cmap = mesh->geometry().cmap();
192 
193  // Get the interpolation points on the reference cells
194  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>& X
195  = element->interpolation_points();
196 
197  mesh->topology_mutable().create_entity_permutations();
198  const std::vector<std::uint32_t>& cell_info
199  = mesh->topology().get_cell_permutation_info();
200 
201  // Push reference coordinates (X) forward to the physical coordinates
202  // (x) for each cell
203  EigenMatrixRowXd x_cell(X.rows(), gdim);
204  std::vector<double> x;
205  EigenMatrixRowXd coordinate_dofs(num_dofs_g, gdim);
206  for (int c = 0; c < num_cells; ++c)
207  {
208  // Get geometry data for current cell
209  auto x_dofs = x_dofmap.links(c);
210  for (int i = 0; i < num_dofs_g; ++i)
211  coordinate_dofs.row(i) = x_g.row(x_dofs[i]).head(gdim);
212 
213  // Push forward coordinates (X -> x)
214  cmap.push_forward(x_cell, X, coordinate_dofs);
215  x.insert(x.end(), x_cell.data(), x_cell.data() + x_cell.size());
216  }
217 
218  // Re-pack points (each row for a given coordinate component) and pad
219  // up to gdim with zero
220  Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor> _x
221  = Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor>::Zero(
222  3, x.size() / gdim);
223  for (int i = 0; i < gdim; ++i)
224  {
225  _x.row(i)
226  = Eigen::Map<Eigen::ArrayXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(
227  x.data() + i, x.size() / gdim,
228  Eigen::InnerStride<Eigen::Dynamic>(gdim));
229  }
230 
231  // Evaluate function at physical points. The returned array has a
232  // number of rows equal to the number of components of the function,
233  // and the number of columns is equal to the number of evaluation
234  // points. Scalar case needs special handling as pybind11 will return
235  // a column array when we need a row array.
236  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values;
237  values = f(_x);
238 
239  // FIXME: This is hack for NumPy/pybind11/Eigen that returns 1D arrays a
240  // column vectors. Fix in the pybind11 layer?
241  if (element->value_size() == 1 and values.rows() > 1)
242  values = values.transpose().eval();
243 
244  // Get dofmap
245  const auto dofmap = u.function_space()->dofmap();
246  assert(dofmap);
247  const int dofmap_bs = dofmap->bs();
248 
249  // NOTE: The below loop over cells could be skipped for some elements,
250  // e.g. Lagrange, where the interpolation is just the identity.
251 
252  // Loop over cells and compute interpolation dofs
253  const int num_scalar_dofs = element->space_dimension() / element_bs;
254  const int value_size = element->value_size() / element_bs;
255 
256  // Check that return type from f is the correct shape
257  if ((values.rows() != value_size * element_bs)
258  || (values.cols() != num_cells * X.rows()))
259  throw std::runtime_error("Interpolation data has the wrong shape.");
260 
261  std::vector<T>& coeffs = u.x()->mutable_array();
262  std::vector<T> _coeffs(num_scalar_dofs);
263  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> _vals;
264  for (int c = 0; c < num_cells; ++c)
265  {
266  auto dofs = dofmap->cell_dofs(c);
267  for (int k = 0; k < element_bs; ++k)
268  {
269  // Extract computed expression values for element block k
270  _vals = values.block(k * value_size, c * X.rows(), value_size, X.rows());
271 
272  // Get element degrees of freedom for block
273  element->interpolate(_vals, cell_info[c], _coeffs);
274  assert(_coeffs.size() == num_scalar_dofs);
275 
276  // Copy interpolation dofs into coefficient vector
277  for (int i = 0; i < num_scalar_dofs; ++i)
278  {
279  const int dof = i * element_bs + k;
280  std::div_t pos = std::div(dof, dofmap_bs);
281  coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
282  }
283  }
284  }
285 }
286 //----------------------------------------------------------------------------
287 template <typename T>
289  Function<T>& u,
290  const std::function<void(
291  Eigen::Ref<
292  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
293  const Eigen::Ref<const Eigen::Array<double, Eigen::Dynamic, 3,
294  Eigen::RowMajor>>&)>& f)
295 {
296  const auto element = u.function_space()->element();
297  assert(element);
298  std::vector<int> vshape(element->value_rank(), 1);
299  for (std::size_t i = 0; i < vshape.size(); ++i)
300  vshape[i] = element->value_dimension(i);
301  const int value_size = std::accumulate(std::begin(vshape), std::end(vshape),
302  1, std::multiplies<>());
303 
304  auto fn =
305  [value_size,
306  &f](const Eigen::Ref<
307  const Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor>>& x) {
308  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values(
309  x.cols(), value_size);
310  f(values, x.transpose());
311  return values;
312  };
313 
314  interpolate<T>(u, fn);
315 }
316 //----------------------------------------------------------------------------
317 
318 } // namespace dolfinx::fem
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:26
void push_forward(Eigen::Ref< Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> x, const Eigen::Ref< const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> &X, const Eigen::Ref< const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> &cell_geometry) const
Compute physical coordinates x for points X in the reference configuration.
Definition: CoordinateElement.cpp:67
This class represents a function in a finite element function space , given by.
Definition: Function.h:43
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:151
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:192
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:46
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:128
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:118
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
void interpolate_c(Function< T > &u, const std::function< void(Eigen::Ref< Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >>, const Eigen::Ref< const Eigen::Array< double, Eigen::Dynamic, 3, Eigen::RowMajor >> &)> &f)
Interpolate an expression f(x)
Definition: interpolate.h:288
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:118