DOLFINx
DOLFINx C++ interface
Function.h
1// Copyright (C) 2003-2022 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 "DofMap.h"
10#include "FiniteElement.h"
11#include "FunctionSpace.h"
12#include "interpolate.h"
13#include <dolfinx/common/IndexMap.h>
14#include <dolfinx/la/Vector.h>
15#include <dolfinx/mesh/Geometry.h>
16#include <dolfinx/mesh/Mesh.h>
17#include <dolfinx/mesh/Topology.h>
18#include <functional>
19#include <memory>
20#include <numeric>
21#include <string>
22#include <utility>
23#include <vector>
24#include <xtensor/xadapt.hpp>
25#include <xtensor/xtensor.hpp>
26#include <xtl/xspan.hpp>
27
28namespace dolfinx::fem
29{
30class FunctionSpace;
31template <typename T>
32class Expression;
33
34template <typename T>
35class Expression;
36
43template <typename T>
45{
46public:
48 using value_type = T;
49
52 explicit Function(std::shared_ptr<const FunctionSpace> V)
53 : _function_space(V),
54 _x(std::make_shared<la::Vector<T>>(V->dofmap()->index_map,
55 V->dofmap()->index_map_bs()))
56 {
57 if (!V->component().empty())
58 {
59 throw std::runtime_error("Cannot create Function from subspace. Consider "
60 "collapsing the function space");
61 }
62 }
63
71 Function(std::shared_ptr<const FunctionSpace> V,
72 std::shared_ptr<la::Vector<T>> x)
73 : _function_space(V), _x(x)
74 {
75 // We do not check for a subspace since this constructor is used for
76 // creating subfunctions
77
78 // Assertion uses '<=' to deal with sub-functions
79 assert(V->dofmap());
80 assert(V->dofmap()->index_map->size_global() * V->dofmap()->index_map_bs()
81 <= _x->bs() * _x->map()->size_global());
82 }
83
84 // Copy constructor
85 Function(const Function& v) = delete;
86
88 Function(Function&& v) = default;
89
91 ~Function() = default;
92
94 Function& operator=(Function&& v) = default;
95
96 // Assignment
97 Function& operator=(const Function& v) = delete;
98
102 Function sub(int i) const
103 {
104 auto sub_space = _function_space->sub({i});
105 assert(sub_space);
106 return Function(sub_space, _x);
107 }
108
113 {
114 // Create new collapsed FunctionSpace
115 auto [V, map] = _function_space->collapse();
116
117 // Create new vector
118 auto x = std::make_shared<la::Vector<T>>(V.dofmap()->index_map,
119 V.dofmap()->index_map_bs());
120
121 // Copy values into new vector
122 xtl::span<const T> x_old = _x->array();
123 xtl::span<T> x_new = x->mutable_array();
124 for (std::size_t i = 0; i < map.size(); ++i)
125 {
126 assert((int)i < x_new.size());
127 assert(map[i] < x_old.size());
128 x_new[i] = x_old[map[i]];
129 }
130
131 return Function(std::make_shared<FunctionSpace>(std::move(V)), x);
132 }
133
136 std::shared_ptr<const FunctionSpace> function_space() const
137 {
138 return _function_space;
139 }
140
142 std::shared_ptr<const la::Vector<T>> x() const { return _x; }
143
145 std::shared_ptr<la::Vector<T>> x() { return _x; }
146
151 const xtl::span<const std::int32_t>& cells)
152 {
153 fem::interpolate(*this, v, cells);
154 }
155
159 {
160 assert(_function_space);
161 assert(_function_space->mesh());
162 int tdim = _function_space->mesh()->topology().dim();
163 auto cell_map = _function_space->mesh()->topology().index_map(tdim);
164 assert(cell_map);
165 std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
166 std::vector<std::int32_t> cells(num_cells, 0);
167 std::iota(cells.begin(), cells.end(), 0);
168
169 fem::interpolate(*this, v, cells);
170 }
171
176 const std::function<xt::xarray<T>(const xt::xtensor<double, 2>&)>& f,
177 const xtl::span<const std::int32_t>& cells)
178 {
179 assert(_function_space);
180 assert(_function_space->element());
181 assert(_function_space->mesh());
182 const xt::xtensor<double, 2> x = fem::interpolation_coords(
183 *_function_space->element(), *_function_space->mesh(), cells);
184 auto fx = f(x);
185 if (int vs = _function_space->element()->value_size();
186 vs == 1 and fx.dimension() == 1)
187 {
188 // Check for scalar-valued functions
189 if (fx.shape(0) != x.shape(1))
190 throw std::runtime_error("Data returned by callable has wrong length");
191 }
192 else
193 {
194 // Check for vector/tensor value
195 if (fx.dimension() != 2)
196 throw std::runtime_error("Expected 2D array of data");
197 if (fx.shape(0) != vs)
198 {
199 throw std::runtime_error(
200 "Data returned by callable has wrong shape(0) size");
201 }
202 if (fx.shape(1) != x.shape(1))
203 {
204 throw std::runtime_error(
205 "Data returned by callable has wrong shape(1) size");
206 }
207 }
208
209 fem::interpolate(*this, fx, cells);
210 }
211
215 const std::function<xt::xarray<T>(const xt::xtensor<double, 2>&)>& f)
216 {
217 assert(_function_space);
218 assert(_function_space->mesh());
219 const int tdim = _function_space->mesh()->topology().dim();
220 auto cell_map = _function_space->mesh()->topology().index_map(tdim);
221 assert(cell_map);
222 std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
223 std::vector<std::int32_t> cells(num_cells, 0);
224 std::iota(cells.begin(), cells.end(), 0);
225 interpolate(f, cells);
226 }
227
235 const xtl::span<const std::int32_t>& cells)
236 {
237 // Check that spaces are compatible
238 assert(_function_space);
239 assert(_function_space->element());
240 std::size_t value_size = e.value_size();
242 throw std::runtime_error("Cannot interpolate Expression with Argument");
243
244 if (value_size != _function_space->element()->value_size())
245 {
246 throw std::runtime_error(
247 "Function value size not equal to Expression value size");
248 }
249
250 if (!xt::allclose(e.X(),
251 _function_space->element()->interpolation_points()))
252 {
253 throw std::runtime_error("Function element interpolation points not "
254 "equal to Expression interpolation points");
255 }
256
257 // Array to hold evaluted Expression
258 std::size_t num_cells = cells.size();
259 std::size_t num_points = e.X().shape(0);
260 xt::xtensor<T, 3> f({num_cells, num_points, value_size});
261
262 // Evaluate Expression at points
263 auto f_view = xt::reshape_view(f, {num_cells, num_points * value_size});
264 e.eval(cells, f_view);
265
266 // Reshape evaluated data to fit interpolate
267 // Expression returns matrix of shape (num_cells, num_points *
268 // value_size), i.e. xyzxyz ordering of dof values per cell per point.
269 // The interpolation uses xxyyzz input, ordered for all points of each
270 // cell, i.e. (value_size, num_cells*num_points)
271 xt::xarray<T> _f = xt::reshape_view(xt::transpose(f, {2, 0, 1}),
272 {value_size, num_cells * num_points});
273
274 // Interpolate values into appropriate space
275 fem::interpolate(*this, _f, cells);
276 }
277
281 {
282 assert(_function_space);
283 assert(_function_space->mesh());
284 const int tdim = _function_space->mesh()->topology().dim();
285 auto cell_map = _function_space->mesh()->topology().index_map(tdim);
286 assert(cell_map);
287 std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
288 std::vector<std::int32_t> cells(num_cells, 0);
289 std::iota(cells.begin(), cells.end(), 0);
290 interpolate(e, cells);
291 }
292
303 void eval(const xt::xtensor<double, 2>& x,
304 const xtl::span<const std::int32_t>& cells,
305 xt::xtensor<T, 2>& u) const
306 {
307 if (cells.empty())
308 return;
309
310 // TODO: This could be easily made more efficient by exploiting points
311 // being ordered by the cell to which they belong.
312
313 if (x.shape(0) != cells.size())
314 {
315 throw std::runtime_error(
316 "Number of points and number of cells must be equal.");
317 }
318 if (x.shape(0) != u.shape(0))
319 {
320 throw std::runtime_error(
321 "Length of array for Function values must be the "
322 "same as the number of points.");
323 }
324
325 // Get mesh
326 assert(_function_space);
327 std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
328 assert(mesh);
329 const std::size_t gdim = mesh->geometry().dim();
330 const std::size_t tdim = mesh->topology().dim();
331 auto map = mesh->topology().index_map(tdim);
332
333 // Get geometry data
335 = mesh->geometry().dofmap();
336 const std::size_t num_dofs_g = mesh->geometry().cmap().dim();
337 xtl::span<const double> x_g = mesh->geometry().x();
338
339 // Get coordinate map
340 const fem::CoordinateElement& cmap = mesh->geometry().cmap();
341
342 // Get element
343 assert(_function_space->element());
344 std::shared_ptr<const fem::FiniteElement> element
345 = _function_space->element();
346 assert(element);
347 const int bs_element = element->block_size();
348 const std::size_t reference_value_size
349 = element->reference_value_size() / bs_element;
350 const std::size_t value_size = element->value_size() / bs_element;
351 const std::size_t space_dimension = element->space_dimension() / bs_element;
352
353 // If the space has sub elements, concatenate the evaluations on the sub
354 // elements
355 const int num_sub_elements = element->num_sub_elements();
356 if (num_sub_elements > 1 and num_sub_elements != bs_element)
357 {
358 throw std::runtime_error("Function::eval is not supported for mixed "
359 "elements. Extract subspaces.");
360 }
361
362 // Create work vector for expansion coefficients
363 std::vector<T> coefficients(space_dimension * bs_element);
364
365 // Get dofmap
366 std::shared_ptr<const fem::DofMap> dofmap = _function_space->dofmap();
367 assert(dofmap);
368 const int bs_dof = dofmap->bs();
369
370 xtl::span<const std::uint32_t> cell_info;
371 if (element->needs_dof_transformations())
372 {
373 mesh->topology_mutable().create_entity_permutations();
374 cell_info = xtl::span(mesh->topology().get_cell_permutation_info());
375 }
376
377 xt::xtensor<double, 2> coordinate_dofs
378 = xt::zeros<double>({num_dofs_g, gdim});
379 xt::xtensor<double, 2> xp = xt::zeros<double>({std::size_t(1), gdim});
380
381 // Loop over points
382 std::fill(u.data(), u.data() + u.size(), 0.0);
383 const xtl::span<const T>& _v = _x->array();
384
385 // -- Lambda function for affine pull-backs
386 xt::xtensor<double, 4> data(cmap.tabulate_shape(1, 1));
387 const xt::xtensor<double, 2> X0(xt::zeros<double>({std::size_t(1), tdim}));
388 cmap.tabulate(1, X0, data);
389 const xt::xtensor<double, 2> dphi_i
390 = xt::view(data, xt::range(1, tdim + 1), 0, xt::all(), 0);
391 auto pull_back_affine = [dphi_i](auto&& X, const auto& cell_geometry,
392 auto&& J, auto&& K, const auto& x) mutable
393 {
394 fem::CoordinateElement::compute_jacobian(dphi_i, cell_geometry, J);
397 X, K, fem::CoordinateElement::x0(cell_geometry), x);
398 };
399
400 xt::xtensor<double, 2> dphi;
401 xt::xtensor<double, 2> X({x.shape(0), tdim});
402 xt::xtensor<double, 3> J = xt::zeros<double>({x.shape(0), gdim, tdim});
403 xt::xtensor<double, 3> K = xt::zeros<double>({x.shape(0), tdim, gdim});
404 xt::xtensor<double, 1> detJ = xt::zeros<double>({x.shape(0)});
405 xt::xtensor<double, 4> phi(cmap.tabulate_shape(1, 1));
406
407 xt::xtensor<double, 2> _Xp({1, tdim});
408 for (std::size_t p = 0; p < cells.size(); ++p)
409 {
410 const int cell_index = cells[p];
411 // Skip negative cell indices
412 if (cell_index < 0)
413 continue;
414
415 // Get cell geometry (coordinate dofs)
416 auto x_dofs = x_dofmap.links(cell_index);
417 assert(x_dofs.size() == num_dofs_g);
418 for (std::size_t i = 0; i < num_dofs_g; ++i)
419 {
420 const int pos = 3 * x_dofs[i];
421 for (std::size_t j = 0; j < gdim; ++j)
422 coordinate_dofs(i, j) = x_g[pos + j];
423 }
424
425 for (std::size_t j = 0; j < gdim; ++j)
426 xp(0, j) = x(p, j);
427
428 auto _J = xt::view(J, p, xt::all(), xt::all());
429 auto _K = xt::view(K, p, xt::all(), xt::all());
430 // Compute reference coordinates X, and J, detJ and K
431 if (cmap.is_affine())
432 {
433 pull_back_affine(_Xp, coordinate_dofs, _J, _K, xp);
434 detJ[p]
436 }
437 else
438 {
439 cmap.pull_back_nonaffine(_Xp, xp, coordinate_dofs);
440 cmap.tabulate(1, _Xp, phi);
441 dphi = xt::view(phi, xt::range(1, tdim + 1), 0, xt::all(), 0);
443 _J);
445 detJ[p]
447 }
448 xt::row(X, p) = xt::row(_Xp, 0);
449 }
450
451 // Prepare basis function data structures
452 xt::xtensor<double, 4> basis_derivatives_reference_values(
453 {1, x.shape(0), space_dimension, reference_value_size});
454 xt::xtensor<double, 3> basis_values(
455 {static_cast<std::size_t>(1), space_dimension, value_size});
456
457 // Compute basis on reference element
458 element->tabulate(basis_derivatives_reference_values, X, 0);
459
460 using u_t = xt::xview<decltype(basis_values)&, std::size_t,
461 xt::xall<std::size_t>, xt::xall<std::size_t>>;
462 using U_t
463 = xt::xview<decltype(basis_derivatives_reference_values)&, std::size_t,
464 std::size_t, xt::xall<std::size_t>, xt::xall<std::size_t>>;
465 using J_t = xt::xview<decltype(J)&, std::size_t, xt::xall<std::size_t>,
466 xt::xall<std::size_t>>;
467 using K_t = xt::xview<decltype(K)&, std::size_t, xt::xall<std::size_t>,
468 xt::xall<std::size_t>>;
469 auto push_forward_fn = element->map_fn<u_t, U_t, J_t, K_t>();
470 const std::function<void(const xtl::span<double>&,
471 const xtl::span<const std::uint32_t>&,
472 std::int32_t, int)>
473 apply_dof_transformation
474 = element->get_dof_transformation_function<double>();
475 const std::size_t num_basis_values = space_dimension * reference_value_size;
476
477 for (std::size_t p = 0; p < cells.size(); ++p)
478 {
479 const int cell_index = cells[p];
480 // Skip negative cell indices
481 if (cell_index < 0)
482 continue;
483
484 // Permute the reference values to account for the cell's orientation
485 apply_dof_transformation(
486 xtl::span(basis_derivatives_reference_values.data()
487 + p * num_basis_values,
488 num_basis_values),
489 cell_info, cell_index, reference_value_size);
490
491 auto _K = xt::view(K, p, xt::all(), xt::all());
492 auto _J = xt::view(J, p, xt::all(), xt::all());
493 auto _U = xt::view(basis_derivatives_reference_values, (std::size_t)0, p,
494 xt::all(), xt::all());
495 auto _u = xt::view(basis_values, (std::size_t)0, xt::all(), xt::all());
496 push_forward_fn(_u, _U, _J, detJ[p], _K);
497
498 // Get degrees of freedom for current cell
499 xtl::span<const std::int32_t> dofs = dofmap->cell_dofs(cell_index);
500 for (std::size_t i = 0; i < dofs.size(); ++i)
501 for (int k = 0; k < bs_dof; ++k)
502 coefficients[bs_dof * i + k] = _v[bs_dof * dofs[i] + k];
503
504 // Compute expansion
505 auto u_row = xt::row(u, p);
506 for (int k = 0; k < bs_element; ++k)
507 {
508 for (std::size_t i = 0; i < space_dimension; ++i)
509 {
510 for (std::size_t j = 0; j < value_size; ++j)
511 {
512 u_row[j * bs_element + k]
513 += coefficients[bs_element * i + k] * basis_values(0, i, j);
514 }
515 }
516 }
517 }
518 }
519
521 std::string name = "u";
522
523private:
524 // The function space
525 std::shared_ptr<const FunctionSpace> _function_space;
526
527 // The vector of expansion coefficients (local)
528 std::shared_ptr<la::Vector<T>> _x;
529};
530} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
Definition: CoordinateElement.h:31
static void compute_jacobian(const U &dphi, const V &cell_geometry, W &&J)
Compute Jacobian for a cell with given geometry using the basis functions and first order derivatives...
Definition: CoordinateElement.h:107
static double compute_jacobian_determinant(const U &J)
Compute the determinant of the Jacobian.
Definition: CoordinateElement.h:130
static void pull_back_affine(xt::xtensor< double, 2 > &X, const xt::xtensor< double, 2 > &K, const std::array< double, 3 > &x0, const xt::xtensor< double, 2 > &x)
Compute reference coordinates X for physical coordinates x for an affine map. For the affine case,...
Definition: CoordinateElement.cpp:89
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition: CoordinateElement.h:116
void pull_back_nonaffine(xt::xtensor< double, 2 > &X, const xt::xtensor< double, 2 > &x, const xt::xtensor< double, 2 > &cell_geometry, double tol=1.0e-8, int maxit=10) const
Compute reference coordinates X for physical coordinates x for a non-affine map.
Definition: CoordinateElement.cpp:106
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Shape of array to fill when calling FiniteElement::tabulate
Definition: CoordinateElement.cpp:44
static std::array< double, 3 > x0(const xt::xtensor< double, 2 > &cell_geometry)
Compute the physical coordinate of the reference point X=(0 , 0, 0)
Definition: CoordinateElement.cpp:81
bool is_affine() const noexcept
Check is geometry map is affine.
Definition: CoordinateElement.h:212
xt::xtensor< double, 4 > tabulate(int nd, const xt::xtensor< double, 2 > &X) const
Evaluate basis values and derivatives at set of points.
Definition: CoordinateElement.cpp:50
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell....
Definition: Expression.h:34
void eval(const xtl::span< const std::int32_t > &cells, U &values) const
Evaluate the expression on cells.
Definition: Expression.h:124
int value_size() const
Get value size.
Definition: Expression.h:215
const xt::xtensor< double, 2 > & X() const
Get evaluation points on reference cell.
Definition: Expression.h:227
std::shared_ptr< const fem::FunctionSpace > argument_function_space() const
Get argument function space.
Definition: Expression.h:81
This class represents a function in a finite element function space , given by.
Definition: Function.h:45
void interpolate(const Expression< T > &e, const xtl::span< const std::int32_t > &cells)
Interpolate an Expression (based on UFL)
Definition: Function.h:234
void interpolate(const Expression< T > &e)
Interpolate an Expression (based on UFL) on all cells.
Definition: Function.h:280
void interpolate(const Function< T > &v, const xtl::span< const std::int32_t > &cells)
Interpolate a Function.
Definition: Function.h:150
void interpolate(const Function< T > &v)
Interpolate a Function.
Definition: Function.h:158
Function collapse() const
Collapse a subfunction (view into a Function) to a stand-alone Function.
Definition: Function.h:112
Function sub(int i) const
Extract subfunction (view into the Function)
Definition: Function.h:102
void eval(const xt::xtensor< double, 2 > &x, const xtl::span< const std::int32_t > &cells, xt::xtensor< T, 2 > &u) const
Evaluate the Function at points.
Definition: Function.h:303
void interpolate(const std::function< xt::xarray< T >(const xt::xtensor< double, 2 > &)> &f)
Interpolate an expression function on the whole domain.
Definition: Function.h:214
std::shared_ptr< const FunctionSpace > function_space() const
Access the function space.
Definition: Function.h:136
std::string name
Name.
Definition: Function.h:521
Function(std::shared_ptr< const FunctionSpace > V)
Create function on given function space.
Definition: Function.h:52
std::shared_ptr< la::Vector< T > > x()
Underlying vector.
Definition: Function.h:145
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:142
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:71
~Function()=default
Destructor.
void interpolate(const std::function< xt::xarray< T >(const xt::xtensor< double, 2 > &)> &f, const xtl::span< const std::int32_t > &cells)
Interpolate an expression function on a list of cells.
Definition: Function.h:175
Function(Function &&v)=default
Move constructor.
Function & operator=(Function &&v)=default
Move assignment.
T value_type
Field type for the Function, e.g. double.
Definition: Function.h:48
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:46
xtl::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:131
Distributed vector.
Definition: Vector.h:25
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
xt::xtensor< double, 2 > interpolation_coords(const fem::FiniteElement &element, const mesh::Mesh &mesh, const xtl::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 xt::xarray< T > &f, const xtl::span< const std::int32_t > &cells)
Interpolate an expression f(x) in a finite element space.
Definition: interpolate.h:383