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