DOLFIN-X
DOLFIN-X C++ interface
Function.h
1 // Copyright (C) 2003-2020 Anders Logg and 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 "interpolate.h"
11 #include <Eigen/Dense>
12 #include <dolfinx/common/IndexMap.h>
13 #include <dolfinx/common/UniqueIdGenerator.h>
14 #include <dolfinx/common/types.h>
15 #include <dolfinx/fem/DofMap.h>
16 #include <dolfinx/fem/FiniteElement.h>
17 #include <dolfinx/la/PETScVector.h>
18 #include <dolfinx/la/Vector.h>
19 #include <dolfinx/mesh/Geometry.h>
20 #include <dolfinx/mesh/Mesh.h>
21 #include <dolfinx/mesh/Topology.h>
22 #include <functional>
23 #include <memory>
24 #include <petscvec.h>
25 #include <string>
26 #include <utility>
27 #include <vector>
28 
29 namespace dolfinx::fem
30 {
31 
32 class FunctionSpace;
33 
40 
41 template <typename T>
42 class Function
43 {
44 public:
47  explicit Function(std::shared_ptr<const FunctionSpace> V)
48  : _id(common::UniqueIdGenerator::id()), _function_space(V),
49  _x(std::make_shared<la::Vector<T>>(V->dofmap()->index_map,
50  V->dofmap()->index_map_bs()))
51  {
52  if (!V->component().empty())
53  {
54  throw std::runtime_error("Cannot create Function from subspace. Consider "
55  "collapsing the function space");
56  }
57  }
58 
65  Function(std::shared_ptr<const FunctionSpace> V,
66  std::shared_ptr<la::Vector<T>> x)
67  : _id(common::UniqueIdGenerator::id()), _function_space(V), _x(x)
68  {
69  // We do not check for a subspace since this constructor is used for
70  // creating subfunctions
71 
72  // Assertion uses '<=' to deal with sub-functions
73  assert(V->dofmap());
74  assert(V->dofmap()->index_map->size_global() * V->dofmap()->index_map_bs()
75  <= _x->bs() * _x->map()->size_global());
76  }
77 
78  // Copy constructor
79  Function(const Function& v) = delete;
80 
83  : name(std::move(v.name)), _id(std::move(v._id)),
84  _function_space(std::move(v._function_space)), _x(std::move(v._x)),
85  _petsc_vector(std::exchange(v._petsc_vector, nullptr))
86  {
87  }
88 
90  virtual ~Function()
91  {
92  if (_petsc_vector)
93  VecDestroy(&_petsc_vector);
94  }
95 
97  Function& operator=(Function&& v) noexcept
98  {
99  name = std::move(v.name);
100  _id = std::move(v._id);
101  _function_space = std::move(v._function_space);
102  _x = std::move(v._x);
103  std::swap(_petsc_vector, v._petsc_vector);
104 
105  return *this;
106  }
107 
108  // Assignment
109  Function& operator=(const Function& v) = delete;
110 
114  Function sub(int i) const
115  {
116  auto sub_space = _function_space->sub({i});
117  assert(sub_space);
118  return Function(sub_space, _x);
119  }
120 
125  {
126  // Create new collapsed FunctionSpace
127  const auto [function_space_new, collapsed_map]
128  = _function_space->collapse();
129 
130  // Create new vector
131  assert(function_space_new);
132  auto vector_new = std::make_shared<la::Vector<T>>(
133  function_space_new->dofmap()->index_map,
134  function_space_new->dofmap()->index_map_bs());
135 
136  // Copy values into new vector
137  const std::vector<T>& x_old = _x->array();
138  std::vector<T>& x_new = vector_new->mutable_array();
139  for (std::size_t i = 0; i < collapsed_map.size(); ++i)
140  {
141  assert((int)i < x_new.size());
142  assert(collapsed_map[i] < x_old.size());
143  x_new[i] = x_old[collapsed_map[i]];
144  }
145 
146  return Function(function_space_new, vector_new);
147  }
148 
151  std::shared_ptr<const FunctionSpace> function_space() const
152  {
153  return _function_space;
154  }
155 
159  Vec vector() const
160  {
161  // Check that this is not a sub function
162  assert(_function_space->dofmap());
163  assert(_function_space->dofmap()->index_map);
164  if (_x->bs() * _x->map()->size_global()
165  != _function_space->dofmap()->index_map->size_global()
166  * _function_space->dofmap()->index_map_bs())
167  {
168  throw std::runtime_error(
169  "Cannot access a non-const vector from a subfunction");
170  }
171 
172  // Check that data type is the same as the PETSc build
173  if constexpr (std::is_same<T, PetscScalar>::value)
174  {
175  if (!_petsc_vector)
176  {
177  _petsc_vector = la::create_ghosted_vector(
178  *_function_space->dofmap()->index_map,
179  _function_space->dofmap()->index_map_bs(), _x->mutable_array());
180  }
181 
182  return _petsc_vector;
183  }
184  else
185  {
186  throw std::runtime_error(
187  "Cannot return PETSc vector wrapper. Type mismatch");
188  }
189  }
190 
192  std::shared_ptr<const la::Vector<T>> x() const { return _x; }
193 
195  std::shared_ptr<la::Vector<T>> x() { return _x; }
196 
199  void interpolate(const Function<T>& v) { fem::interpolate(*this, v); }
200 
203  void
204  interpolate(const std::function<
205  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
206  const Eigen::Ref<const Eigen::Array<double, 3, Eigen::Dynamic,
207  Eigen::RowMajor>>&)>& f)
208  {
209  fem::interpolate(*this, f);
210  }
211 
222  void
223  eval(const Eigen::Ref<
224  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>>& x,
225  const tcb::span<const std::int32_t>& cells,
226  Eigen::Ref<
227  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
228  u) const
229  {
230  // TODO: This could be easily made more efficient by exploiting points
231  // being ordered by the cell to which they belong.
232 
233  if (x.rows() != (int)cells.size())
234  {
235  throw std::runtime_error(
236  "Number of points and number of cells must be equal.");
237  }
238  if (x.rows() != u.rows())
239  {
240  throw std::runtime_error(
241  "Length of array for Function values must be the "
242  "same as the number of points.");
243  }
244 
245  // Get mesh
246  assert(_function_space);
247  std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
248  assert(mesh);
249  const int gdim = mesh->geometry().dim();
250  const int tdim = mesh->topology().dim();
251  auto map = mesh->topology().index_map(tdim);
252 
253  // Get geometry data
254  const graph::AdjacencyList<std::int32_t>& x_dofmap
255  = mesh->geometry().dofmap();
256  // FIXME: Add proper interface for num coordinate dofs
257  const int num_dofs_g = x_dofmap.num_links(0);
258  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
259  = mesh->geometry().x();
260 
261  // Get coordinate map
262  const fem::CoordinateElement& cmap = mesh->geometry().cmap();
263 
264  // Get element
265  assert(_function_space->element());
266  std::shared_ptr<const fem::FiniteElement> element
267  = _function_space->element();
268  assert(element);
269  const int bs_element = element->block_size();
270  const int reference_value_size
271  = element->reference_value_size() / bs_element;
272  const int value_size = element->value_size() / bs_element;
273  const int space_dimension = element->space_dimension() / bs_element;
274 
275  // If the space has sub elements, concatenate the evaluations on the sub
276  // elements
277  const int num_sub_elements = element->num_sub_elements();
278  if (num_sub_elements > 1 and num_sub_elements != bs_element)
279  {
280  throw std::runtime_error("Function::eval is not supported for mixed "
281  "elements. Extract subspaces.");
282  }
283 
284  // Prepare geometry data structures
285  Eigen::Tensor<double, 3, Eigen::RowMajor> J(1, gdim, tdim);
286  std::array<double, 1> detJ;
287  Eigen::Tensor<double, 3, Eigen::RowMajor> K(1, tdim, gdim);
288  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> X(
289  1, tdim);
290 
291  // Prepare basis function data structures
292  Eigen::Tensor<double, 3, Eigen::RowMajor> basis_reference_values(
293  1, space_dimension, reference_value_size);
294  Eigen::Tensor<double, 3, Eigen::RowMajor> basis_values(1, space_dimension,
295  value_size);
296 
297  // Create work vector for expansion coefficients
298  Eigen::Matrix<T, 1, Eigen::Dynamic> coefficients(space_dimension
299  * bs_element);
300 
301  // Get dofmap
302  std::shared_ptr<const fem::DofMap> dofmap = _function_space->dofmap();
303  assert(dofmap);
304  const int bs_dof = dofmap->bs();
305 
306  mesh->topology_mutable().create_entity_permutations();
307  const std::vector<std::uint32_t>& cell_info
308  = mesh->topology().get_cell_permutation_info();
309  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
310  coordinate_dofs(num_dofs_g, gdim);
311 
312  // Loop over points
313  u.setZero();
314  const std::vector<T>& _v = _x->mutable_array();
315  for (std::size_t p = 0; p < cells.size(); ++p)
316  {
317  const int cell_index = cells[p];
318 
319  // Skip negative cell indices
320  if (cell_index < 0)
321  continue;
322 
323  // Get cell geometry (coordinate dofs)
324  auto x_dofs = x_dofmap.links(cell_index);
325  for (int i = 0; i < num_dofs_g; ++i)
326  coordinate_dofs.row(i) = x_g.row(x_dofs[i]).head(gdim);
327 
328  // Compute reference coordinates X, and J, detJ and K
329  cmap.compute_reference_geometry(X, J, detJ, K, x.row(p).head(gdim),
330  coordinate_dofs);
331 
332  // Compute basis on reference element
333  element->evaluate_reference_basis(basis_reference_values, X);
334 
335  // Permute the reference values to account for the cell's orientation
336  element->apply_dof_transformation(basis_reference_values.data(),
337  cell_info[cell_index],
338  reference_value_size);
339 
340  // Push basis forward to physical element
341  element->transform_reference_basis(basis_values, basis_reference_values,
342  X, J, detJ, K);
343 
344  // Get degrees of freedom for current cell
345  tcb::span<const std::int32_t> dofs = dofmap->cell_dofs(cell_index);
346  for (std::size_t i = 0; i < dofs.size(); ++i)
347  for (int k = 0; k < bs_dof; ++k)
348  coefficients[bs_dof * i + k] = _v[bs_dof * dofs[i] + k];
349 
350  // Compute expansion
351  auto u_row = u.row(p);
352  for (int k = 0; k < bs_element; ++k)
353  {
354  for (int i = 0; i < space_dimension; ++i)
355  {
356  for (int j = 0; j < value_size; ++j)
357  {
358  // TODO: Find an Eigen shortcut for this operation?
359  u_row[j * bs_element + k]
360  += coefficients[bs_element * i + k] * basis_values(0, i, j);
361  }
362  }
363  }
364  }
365  }
366 
369  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
371  {
372  assert(_function_space);
373  std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
374  assert(mesh);
375  const int tdim = mesh->topology().dim();
376 
377  // Compute in tensor (one for scalar function, . . .)
378  const int value_size_loc = _function_space->element()->value_size();
379 
380  // Resize Array for holding point values
381  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
382  point_values(mesh->geometry().x().rows(), value_size_loc);
383 
384  // Prepare cell geometry
385  const graph::AdjacencyList<std::int32_t>& x_dofmap
386  = mesh->geometry().dofmap();
387 
388  // FIXME: Add proper interface for num coordinate dofs
389  const int num_dofs_g = x_dofmap.num_links(0);
390  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
391  = mesh->geometry().x();
392 
393  // Interpolate point values on each cell (using last computed value if
394  // not continuous, e.g. discontinuous Galerkin methods)
395  auto map = mesh->topology().index_map(tdim);
396  assert(map);
397  const std::int32_t num_cells = map->size_local() + map->num_ghosts();
398 
399  std::vector<std::int32_t> cells(x_g.rows());
400  for (std::int32_t c = 0; c < num_cells; ++c)
401  {
402  // Get coordinates for all points in cell
403  tcb::span<const std::int32_t> dofs = x_dofmap.links(c);
404  for (int i = 0; i < num_dofs_g; ++i)
405  cells[dofs[i]] = c;
406  }
407 
408  eval(x_g, cells, point_values);
409 
410  return point_values;
411  }
412 
414  std::string name = "u";
415 
417  std::size_t id() const { return _id; }
418 
419 private:
420  // ID
421  std::size_t _id;
422 
423  // The function space
424  std::shared_ptr<const FunctionSpace> _function_space;
425 
426  // The vector of expansion coefficients (local)
427  std::shared_ptr<la::Vector<T>> _x;
428 
429  // PETSc wrapper of the expansion coefficients
430  mutable Vec _petsc_vector = nullptr;
431 };
432 } // namespace dolfinx::fem
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:26
void compute_reference_geometry(Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &X, Eigen::Tensor< double, 3, Eigen::RowMajor > &J, tcb::span< double > detJ, Eigen::Tensor< double, 3, Eigen::RowMajor > &K, 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 reference coordinates X, and J, detJ and K for physical coordinates x.
Definition: CoordinateElement.cpp:85
This class represents a function in a finite element function space , given by.
Definition: Function.h:43
Vec vector() const
Return vector of expansion coefficients as a PETSc Vec. Throws an exception a PETSc Vec cannot be cre...
Definition: Function.h:159
Function(std::shared_ptr< const FunctionSpace > V, std::shared_ptr< la::Vector< T >> x)
Create function on given function space with a given vector.
Definition: Function.h:65
void interpolate(const Function< T > &v)
Interpolate a Function (on possibly non-matching meshes)
Definition: Function.h:199
Function & operator=(Function &&v) noexcept
Move assignment.
Definition: Function.h:97
Function(Function &&v)
Move constructor.
Definition: Function.h:82
void eval(const Eigen::Ref< const Eigen::Array< double, Eigen::Dynamic, 3, Eigen::RowMajor >> &x, const tcb::span< const std::int32_t > &cells, Eigen::Ref< Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> u) const
Evaluate the Function at points.
Definition: Function.h:223
std::shared_ptr< la::Vector< T > > x()
Underlying vector.
Definition: Function.h:195
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:151
virtual ~Function()
Destructor.
Definition: Function.h:90
void interpolate(const std::function< Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >(const Eigen::Ref< const Eigen::Array< double, 3, Eigen::Dynamic, Eigen::RowMajor >> &)> &f)
Interpolate an expression.
Definition: Function.h:204
std::size_t id() const
ID.
Definition: Function.h:417
std::string name
Name.
Definition: Function.h:414
Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > compute_point_values() const
Compute values at all mesh 'nodes'.
Definition: Function.h:370
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:192
Function collapse() const
Collapse a subfunction (view into the Function) to a stand-alone Function.
Definition: Function.h:124
Function(std::shared_ptr< const FunctionSpace > V)
Create function on given function space.
Definition: Function.h:47
Function sub(int i) const
Extract subfunction (view into the Function)
Definition: Function.h:114
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
Distributed vector.
Definition: Vector.h:19
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
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
Vec create_ghosted_vector(const common::IndexMap &map, int bs, tcb::span< PetscScalar > x)
Create a PETSc Vec that wraps the data in an array.
Definition: PETScVector.cpp:65