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