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::function
18 {
19 
20 template <typename T>
21 class Function;
22 
26 template <typename T>
27 void interpolate(Function<T>& u, const Function<T>& v);
28 
32 template <typename T>
33 void interpolate(
34  Function<T>& u,
35  const std::function<
36  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
37  const Eigen::Ref<const Eigen::Array<double, 3, Eigen::Dynamic,
38  Eigen::RowMajor>>&)>& f);
39 
48 template <typename T>
49 void interpolate_c(
50  Function<T>& u,
51  const std::function<void(
52  Eigen::Ref<
53  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
54  const Eigen::Ref<const Eigen::Array<double, Eigen::Dynamic, 3,
55  Eigen::RowMajor>>&)>& f);
56 
57 namespace detail
58 {
59 
60 // Interpolate data. Fills coefficients using 'values', which are the
61 // values of an expression at each dof.
62 template <typename T>
63 void interpolate_values(
64  Function<T>& u,
65  const Eigen::Ref<const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic,
66  Eigen::RowMajor>>& values)
67 {
68  Eigen::Matrix<T, Eigen::Dynamic, 1>& coefficients = u.x()->array();
69  coefficients = Eigen::Map<const Eigen::Array<T, Eigen::Dynamic, 1>>(
70  values.data(), coefficients.rows());
71 }
72 
73 template <typename T>
74 void interpolate_from_any(Function<T>& u, const Function<T>& v)
75 {
76  assert(v.function_space());
77  const auto element = u.function_space()->element();
78  assert(element);
79  if (!v.function_space()->has_element(*element))
80  {
81  throw std::runtime_error("Restricting finite elements function in "
82  "different elements not supported.");
83  }
84 
85  const auto mesh = u.function_space()->mesh();
86  assert(mesh);
87  assert(v.function_space()->mesh());
88  if (mesh->id() != v.function_space()->mesh()->id())
89  {
90  throw std::runtime_error(
91  "Interpolation on different meshes not supported (yet).");
92  }
93  const int tdim = mesh->topology().dim();
94 
95  // Get dofmaps
96  assert(v.function_space());
97  std::shared_ptr<const fem::DofMap> dofmap_v = v.function_space()->dofmap();
98  assert(dofmap_v);
99  auto map = mesh->topology().index_map(tdim);
100  assert(map);
101 
102  Eigen::Matrix<T, Eigen::Dynamic, 1>& expansion_coefficients = u.x()->array();
103 
104  // Iterate over mesh and interpolate on each cell
105  const auto dofmap_u = u.function_space()->dofmap();
106  const Eigen::Matrix<T, Eigen::Dynamic, 1>& v_array = v.x()->array();
107  const int num_cells = map->size_local() + map->num_ghosts();
108  for (int c = 0; c < num_cells; ++c)
109  {
110  auto dofs_v = dofmap_v->cell_dofs(c);
111  auto cell_dofs = dofmap_u->cell_dofs(c);
112  assert(dofs_v.size() == cell_dofs.size());
113  for (Eigen::Index i = 0; i < dofs_v.size(); ++i)
114  expansion_coefficients[cell_dofs[i]] = v_array[dofs_v[i]];
115  }
116 }
117 
118 } // namespace detail
119 
120 //----------------------------------------------------------------------------
121 template <typename T>
123 {
124  assert(u.function_space());
125  const auto element = u.function_space()->element();
126  assert(element);
127 
128  // Check that function ranks match
129  if (int rank_v = v.function_space()->element()->value_rank();
130  element->value_rank() != rank_v)
131  {
132  throw std::runtime_error("Cannot interpolate function into function space. "
133  "Rank of function ("
134  + std::to_string(rank_v)
135  + ") does not match rank of function space ("
136  + std::to_string(element->value_rank()) + ")");
137  }
138 
139  // Check that function dimension match
140  for (int i = 0; i < element->value_rank(); ++i)
141  {
142  if (int v_dim = v.function_space()->element()->value_dimension(i);
143  element->value_dimension(i) != v_dim)
144  {
145  throw std::runtime_error(
146  "Cannot interpolate function into function space. "
147  "Dimension "
148  + std::to_string(i) + " of function (" + std::to_string(v_dim)
149  + ") does not match dimension " + std::to_string(i)
150  + " of function space(" + std::to_string(element->value_dimension(i))
151  + ")");
152  }
153  }
154 
155  detail::interpolate_from_any(u, v);
156 }
157 //----------------------------------------------------------------------------
158 template <typename T>
160  Function<T>& u,
161  const std::function<
162  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
163  const Eigen::Ref<const Eigen::Array<double, 3, Eigen::Dynamic,
164  Eigen::RowMajor>>&)>& f)
165 {
166  // u.function_space()->interpolate(u.x()->array(), f);
167  assert(u.function_space());
168 
169  // Evaluate expression at dof points
170  const Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor> x
171  = u.function_space()
172  ->tabulate_scalar_subspace_dof_coordinates()
173  .transpose();
174  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values
175  = f(x);
176 
177  const auto element = u.function_space()->element();
178  assert(element);
179  std::vector<int> vshape(element->value_rank(), 1);
180  for (std::size_t i = 0; i < vshape.size(); ++i)
181  vshape[i] = element->value_dimension(i);
182  const int value_size = std::accumulate(std::begin(vshape), std::end(vshape),
183  1, std::multiplies<>());
184 
185  auto mesh = u.function_space()->mesh();
186  assert(mesh);
187  const int tdim = mesh->topology().dim();
188  if (element->family() == "Mixed")
189  {
190  // Extract the correct components of the result for each subelement of the
191  // MixedElement
192  std::shared_ptr<const fem::DofMap> dofmap = u.function_space()->dofmap();
193  auto map = mesh->topology().index_map(tdim);
194  assert(map);
195  const int num_cells = map->size_local() + map->num_ghosts();
196 
197  Eigen::Array<T, 1, Eigen::Dynamic> mixed_values(values.cols());
198  int value_offset = 0;
199  for (int i = 0; i < element->num_sub_elements(); ++i)
200  {
201  const std::vector<int> component = {i};
202  const fem::DofMap sub_dofmap = dofmap->extract_sub_dofmap(component);
203  std::shared_ptr<const fem::FiniteElement> sub_element
204  = element->extract_sub_element(component);
205  const int element_block_size = sub_element->block_size();
206  for (std::size_t cell = 0; cell < num_cells; ++cell)
207  {
208  const auto cell_dofs = sub_dofmap.cell_dofs(cell);
209  const std::size_t scalar_dofs = cell_dofs.rows() / element_block_size;
210  for (std::size_t dof = 0; dof < scalar_dofs; ++dof)
211  for (int b = 0; b < element_block_size; ++b)
212  mixed_values(cell_dofs[element_block_size * dof + b])
213  = values(value_offset + b, cell_dofs[element_block_size * dof]);
214  }
215  value_offset += element_block_size;
216  }
217  detail::interpolate_values<T>(u, mixed_values);
218  return;
219  }
220 
221  // Note: pybind11 maps 1D NumPy arrays to column vectors for
222  // Eigen::Array<T, Eigen::Dynamic,Eigen::Dynamic, Eigen::RowMajor>
223  // types, therefore we need to handle vectors as a special case.
224  if (values.cols() == 1 and values.rows() != 1)
225  {
226  if (values.rows() != x.cols())
227  {
228  throw std::runtime_error("Number of computed values is not equal to the "
229  "number of evaluation points. (1)");
230  }
231  detail::interpolate_values<T>(u, values);
232  }
233  else
234  {
235  if (values.rows() != value_size)
236  throw std::runtime_error("Values shape is incorrect. (2)");
237  if (values.cols() != x.cols())
238  {
239  throw std::runtime_error("Number of computed values is not equal to the "
240  "number of evaluation points. (2)");
241  }
242 
243  detail::interpolate_values<T>(u, values.transpose());
244  }
245 }
246 //----------------------------------------------------------------------------
247 template <typename T>
249  Function<T>& u,
250  const std::function<void(
251  Eigen::Ref<
252  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
253  const Eigen::Ref<const Eigen::Array<double, Eigen::Dynamic, 3,
254  Eigen::RowMajor>>&)>& f)
255 {
256  // Build list of points at which to evaluate the Expression
257  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor> x
258  = u.function_space()->tabulate_dof_coordinates();
259 
260  // Evaluate expression at points
261  const auto element = u.function_space()->element();
262  assert(element);
263  std::vector<int> vshape(element->value_rank(), 1);
264  for (std::size_t i = 0; i < vshape.size(); ++i)
265  vshape[i] = element->value_dimension(i);
266  const int value_size = std::accumulate(std::begin(vshape), std::end(vshape),
267  1, std::multiplies<>());
268  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values(
269  x.rows(), value_size);
270  f(values, x);
271 
272  detail::interpolate_values<T>(u, values);
273 }
274 //----------------------------------------------------------------------------
275 
276 } // namespace dolfinx::function
Degree-of-freedom map.
Definition: DofMap.h:66
Eigen::Array< std::int32_t, Eigen::Dynamic, 1 >::ConstSegmentReturnType cell_dofs(int cell) const
Local-to-global mapping of dofs on a cell.
Definition: DofMap.h:95
DofMap extract_sub_dofmap(const std::vector< int > &component) const
Extract subdofmap component.
Definition: DofMap.cpp:207
This class represents a function in a finite element function space , given by.
Definition: Function.h:42
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:136
Functions tools, including FEM functions and pointwise defined functions.
Definition: assembler.h:19
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). This interface uses an expression function f that has an in/out argum...
Definition: interpolate.h:248
void interpolate(Function< T > &u, const Function< T > &v)
Interpolate a Function (on possibly non-matching meshes)
Definition: interpolate.h:122