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 <span>
22#include <string>
23#include <utility>
24#include <vector>
25#include <xtensor/xadapt.hpp>
26#include <xtensor/xtensor.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 std::span<const T> x_old = _x->array();
123 std::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
150 void interpolate(const Function<T>& v, std::span<const std::int32_t> cells)
151 {
152 fem::interpolate(*this, v, cells);
153 }
154
158 {
159 assert(_function_space);
160 assert(_function_space->mesh());
161 int tdim = _function_space->mesh()->topology().dim();
162 auto cell_map = _function_space->mesh()->topology().index_map(tdim);
163 assert(cell_map);
164 std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
165 std::vector<std::int32_t> cells(num_cells, 0);
166 std::iota(cells.begin(), cells.end(), 0);
167
168 fem::interpolate(*this, v, cells);
169 }
170
175 const std::function<xt::xarray<T>(const xt::xtensor<double, 2>&)>& f,
176 std::span<const std::int32_t> cells)
177 {
178 assert(_function_space);
179 assert(_function_space->element());
180 assert(_function_space->mesh());
181 const std::vector<double> x = fem::interpolation_coords(
182 *_function_space->element(), *_function_space->mesh(), cells);
183 auto _x = xt::adapt(x, std::vector<std::size_t>{3, x.size() / 3});
184 auto fx = f(_x);
185 assert(fx.dimension() <= 2);
186 if (int vs = _function_space->element()->value_size();
187 vs == 1 and fx.dimension() == 1)
188 {
189 // Check for scalar-valued functions
190 if (fx.shape(0) != x.size() / 3)
191 throw std::runtime_error("Data returned by callable has wrong length");
192 }
193 else
194 {
195 // Check for vector/tensor value
196 if (fx.dimension() != 2)
197 throw std::runtime_error("Expected 2D array of data");
198 if (fx.shape(0) != vs)
199 {
200 throw std::runtime_error(
201 "Data returned by callable has wrong shape(0) size");
202 }
203 if (fx.shape(1) != x.size() / 3)
204 {
205 throw std::runtime_error(
206 "Data returned by callable has wrong shape(1) size");
207 }
208 }
209
210 std::array<std::size_t, 2> fshape;
211 if (fx.dimension() == 1)
212 fshape = {1, fx.shape(0)};
213 else
214 fshape = {fx.shape(0), fx.shape(1)};
215
216 fem::interpolate(*this, std::span<const T>(fx.data(), fx.size()), fshape,
217 cells);
218 }
219
223 const std::function<xt::xarray<T>(const xt::xtensor<double, 2>&)>& f)
224 {
225 assert(_function_space);
226 assert(_function_space->mesh());
227 const int tdim = _function_space->mesh()->topology().dim();
228 auto cell_map = _function_space->mesh()->topology().index_map(tdim);
229 assert(cell_map);
230 std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
231 std::vector<std::int32_t> cells(num_cells, 0);
232 std::iota(cells.begin(), cells.end(), 0);
233 interpolate(f, cells);
234 }
235
242 void interpolate(const Expression<T>& e, std::span<const std::int32_t> cells)
243 {
244 // Check that spaces are compatible
245 assert(_function_space);
246 assert(_function_space->element());
247 std::size_t value_size = e.value_size();
249 throw std::runtime_error("Cannot interpolate Expression with Argument");
250
251 if (value_size != _function_space->element()->value_size())
252 {
253 throw std::runtime_error(
254 "Function value size not equal to Expression value size");
255 }
256
257 {
258 // Compatibility check
259 auto [X0, shape0] = e.X();
260 auto [X1, shape1] = _function_space->element()->interpolation_points();
261 if (shape0 != shape1)
262 {
263 throw std::runtime_error(
264 "Function element interpolation points has different shape to "
265 "Expression interpolation points");
266 }
267 for (std::size_t i = 0; i < X0.size(); ++i)
268 {
269 if (std::abs(X0[i] - X1[i]) > 1.0e-10)
270 {
271 throw std::runtime_error("Function element interpolation points not "
272 "equal to Expression interpolation points");
273 }
274 }
275 }
276
277 // Array to hold evaluted Expression
278 std::size_t num_cells = cells.size();
279 std::size_t num_points = e.X().second[0];
280 xt::xtensor<T, 3> f({num_cells, num_points, value_size});
281
282 // Evaluate Expression at points
283 auto f_view = xt::reshape_view(f, {num_cells, num_points * value_size});
284 e.eval(cells, f_view);
285
286 // Reshape evaluated data to fit interpolate
287 // Expression returns matrix of shape (num_cells, num_points *
288 // value_size), i.e. xyzxyz ordering of dof values per cell per point.
289 // The interpolation uses xxyyzz input, ordered for all points of each
290 // cell, i.e. (value_size, num_cells*num_points)
291 xt::xarray<T> _f = xt::reshape_view(xt::transpose(f, {2, 0, 1}),
292 {value_size, num_cells * num_points});
293
294 // Interpolate values into appropriate space
295 fem::interpolate(*this, std::span<const T>(_f.data(), _f.size()),
296 {_f.shape(0), _f.shape(1)}, cells);
297 }
298
302 {
303 assert(_function_space);
304 assert(_function_space->mesh());
305 const int tdim = _function_space->mesh()->topology().dim();
306 auto cell_map = _function_space->mesh()->topology().index_map(tdim);
307 assert(cell_map);
308 std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
309 std::vector<std::int32_t> cells(num_cells, 0);
310 std::iota(cells.begin(), cells.end(), 0);
311 interpolate(e, cells);
312 }
313
326 void eval(std::span<const double> x, std::array<std::size_t, 2> xshape,
327 std::span<const std::int32_t> cells, std::span<T> u,
328 std::array<std::size_t, 2> ushape) const
329 {
330 if (cells.empty())
331 return;
332
333 assert(x.size() == xshape[0] * xshape[1]);
334 assert(u.size() == ushape[0] * ushape[1]);
335
336 // TODO: This could be easily made more efficient by exploiting points
337 // being ordered by the cell to which they belong.
338
339 if (xshape[0] != cells.size())
340 {
341 throw std::runtime_error(
342 "Number of points and number of cells must be equal.");
343 }
344 if (xshape[0] != ushape[0])
345 {
346 throw std::runtime_error(
347 "Length of array for Function values must be the "
348 "same as the number of points.");
349 }
350
351 // Get mesh
352 assert(_function_space);
353 std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
354 assert(mesh);
355 const std::size_t gdim = mesh->geometry().dim();
356 const std::size_t tdim = mesh->topology().dim();
357 auto map = mesh->topology().index_map(tdim);
358
359 // Get geometry data
361 = mesh->geometry().dofmap();
362 const std::size_t num_dofs_g = mesh->geometry().cmap().dim();
363 std::span<const double> x_g = mesh->geometry().x();
364
365 // Get coordinate map
366 const CoordinateElement& cmap = mesh->geometry().cmap();
367
368 // Get element
369 std::shared_ptr<const FiniteElement> element = _function_space->element();
370 assert(element);
371 const int bs_element = element->block_size();
372 const std::size_t reference_value_size
373 = element->reference_value_size() / bs_element;
374 const std::size_t value_size = element->value_size() / bs_element;
375 const std::size_t space_dimension = element->space_dimension() / bs_element;
376
377 // If the space has sub elements, concatenate the evaluations on the
378 // sub elements
379 const int num_sub_elements = element->num_sub_elements();
380 if (num_sub_elements > 1 and num_sub_elements != bs_element)
381 {
382 throw std::runtime_error("Function::eval is not supported for mixed "
383 "elements. Extract subspaces.");
384 }
385
386 // Create work vector for expansion coefficients
387 std::vector<T> coefficients(space_dimension * bs_element);
388
389 // Get dofmap
390 std::shared_ptr<const DofMap> dofmap = _function_space->dofmap();
391 assert(dofmap);
392 const int bs_dof = dofmap->bs();
393
394 std::span<const std::uint32_t> cell_info;
395 if (element->needs_dof_transformations())
396 {
397 mesh->topology_mutable().create_entity_permutations();
398 cell_info = std::span(mesh->topology().get_cell_permutation_info());
399 }
400
401 namespace stdex = std::experimental;
402 using cmdspan4_t
403 = stdex::mdspan<const double, stdex::dextents<std::size_t, 4>>;
404 using mdspan2_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
405 using mdspan3_t = stdex::mdspan<double, stdex::dextents<std::size_t, 3>>;
406
407 std::vector<double> coord_dofs_b(num_dofs_g * gdim);
408 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
409 std::vector<double> xp_b(1 * gdim);
410 mdspan2_t xp(xp_b.data(), 1, gdim);
411
412 // Loop over points
413 std::fill(u.data(), u.data() + u.size(), 0.0);
414 const std::span<const T>& _v = _x->array();
415
416 // Evaluate geometry basis at point (0, 0, 0) on the reference cell.
417 // Used in affine case.
418 std::array<std::size_t, 4> phi0_shape = cmap.tabulate_shape(1, 1);
419 std::vector<double> phi0_b(std::reduce(phi0_shape.begin(), phi0_shape.end(),
420 1, std::multiplies{}));
421 cmdspan4_t phi0(phi0_b.data(), phi0_shape);
422 cmap.tabulate(1, std::vector<double>(tdim), {1, tdim}, phi0_b);
423 auto dphi0 = stdex::submdspan(phi0, std::pair(1, tdim + 1), 0,
424 stdex::full_extent, 0);
425
426 // Data structure for evaluating geometry basis at specific points.
427 // Used in non-affine case.
428 std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(1, 1);
429 std::vector<double> phi_b(
430 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
431 cmdspan4_t phi(phi_b.data(), phi_shape);
432 auto dphi = stdex::submdspan(phi, std::pair(1, tdim + 1), 0,
433 stdex::full_extent, 0);
434
435 // Reference coordinates for each point
436 std::vector<double> Xb(xshape[0] * tdim);
437 mdspan2_t X(Xb.data(), xshape[0], tdim);
438
439 // Geometry data at each point
440 std::vector<double> J_b(xshape[0] * gdim * tdim);
441 mdspan3_t J(J_b.data(), xshape[0], gdim, tdim);
442 std::vector<double> K_b(xshape[0] * tdim * gdim);
443 mdspan3_t K(K_b.data(), xshape[0], tdim, gdim);
444 std::vector<double> detJ(xshape[0]);
445 std::vector<double> det_scratch(2 * gdim * tdim);
446
447 // Prepare geometry data in each cell
448 for (std::size_t p = 0; p < cells.size(); ++p)
449 {
450 const int cell_index = cells[p];
451
452 // Skip negative cell indices
453 if (cell_index < 0)
454 continue;
455
456 // Get cell geometry (coordinate dofs)
457 auto x_dofs = x_dofmap.links(cell_index);
458 assert(x_dofs.size() == num_dofs_g);
459 for (std::size_t i = 0; i < num_dofs_g; ++i)
460 {
461 const int pos = 3 * x_dofs[i];
462 for (std::size_t j = 0; j < gdim; ++j)
463 coord_dofs(i, j) = x_g[pos + j];
464 }
465
466 for (std::size_t j = 0; j < gdim; ++j)
467 xp(0, j) = x[p * xshape[1] + j];
468
469 auto _J = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
470 auto _K = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
471
472 std::array<double, 3> Xpb = {0, 0, 0};
473 stdex::mdspan<double,
474 stdex::extents<std::size_t, 1, stdex::dynamic_extent>>
475 Xp(Xpb.data(), 1, tdim);
476
477 // Compute reference coordinates X, and J, detJ and K
478 if (cmap.is_affine())
479 {
480 CoordinateElement::compute_jacobian(dphi0, coord_dofs, _J);
482 std::array<double, 3> x0 = {0, 0, 0};
483 for (std::size_t i = 0; i < coord_dofs.extent(1); ++i)
484 x0[i] += coord_dofs(0, i);
486 detJ[p]
488 }
489 else
490 {
491 // Pull-back physical point xp to reference coordinate Xp
492 cmap.pull_back_nonaffine(Xp, xp, coord_dofs);
493
494 cmap.tabulate(1, std::span(Xpb.data(), tdim), {1, tdim}, phi_b);
495 CoordinateElement::compute_jacobian(dphi, coord_dofs, _J);
497 detJ[p]
499 }
500
501 for (std::size_t j = 0; j < X.extent(1); ++j)
502 X(p, j) = Xpb[j];
503 }
504
505 // Prepare basis function data structures
506 std::vector<double> basis_derivatives_reference_values_b(
507 1 * xshape[0] * space_dimension * reference_value_size);
508 cmdspan4_t basis_derivatives_reference_values(
509 basis_derivatives_reference_values_b.data(), 1, xshape[0],
510 space_dimension, reference_value_size);
511 std::vector<double> basis_values_b(space_dimension * value_size);
512 mdspan2_t basis_values(basis_values_b.data(), space_dimension, value_size);
513
514 // Compute basis on reference element
515 element->tabulate(basis_derivatives_reference_values_b, Xb,
516 {X.extent(0), X.extent(1)}, 0);
517
518 using xu_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
519 using xU_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
520 using xJ_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
521 using xK_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
522 auto push_forward_fn
523 = element->basix_element().map_fn<xu_t, xU_t, xJ_t, xK_t>();
524
525 auto apply_dof_transformation
526 = element->get_dof_transformation_function<double>();
527 const std::size_t num_basis_values = space_dimension * reference_value_size;
528
529 for (std::size_t p = 0; p < cells.size(); ++p)
530 {
531 const int cell_index = cells[p];
532
533 // Skip negative cell indices
534 if (cell_index < 0)
535 continue;
536
537 // Permute the reference values to account for the cell's
538 // orientation
539 apply_dof_transformation(
540 std::span(basis_derivatives_reference_values_b.data()
541 + p * num_basis_values,
542 num_basis_values),
543 cell_info, cell_index, reference_value_size);
544
545 {
546 auto _U = stdex::submdspan(basis_derivatives_reference_values, 0, p,
547 stdex::full_extent, stdex::full_extent);
548 auto _J
549 = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
550 auto _K
551 = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
552 push_forward_fn(basis_values, _U, _J, detJ[p], _K);
553 }
554
555 // Get degrees of freedom for current cell
556 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell_index);
557 for (std::size_t i = 0; i < dofs.size(); ++i)
558 for (int k = 0; k < bs_dof; ++k)
559 coefficients[bs_dof * i + k] = _v[bs_dof * dofs[i] + k];
560
561 // Compute expansion
562 for (int k = 0; k < bs_element; ++k)
563 {
564 for (std::size_t i = 0; i < space_dimension; ++i)
565 {
566 for (std::size_t j = 0; j < value_size; ++j)
567 {
568 u[p * ushape[1] + (j * bs_element + k)]
569 += coefficients[bs_element * i + k] * basis_values(i, j);
570 }
571 }
572 }
573 }
574 }
575
577 std::string name = "u";
578
579private:
580 // The function space
581 std::shared_ptr<const FunctionSpace> _function_space;
582
583 // The vector of expansion coefficients (local)
584 std::shared_ptr<la::Vector<T>> _x;
585};
586} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
A CoordinateElement manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:32
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:100
static void pull_back_affine(U &&X, const V &K, const std::array< double, 3 > &x0, const W &x)
Compute reference coordinates X for physical coordinates x for an affine map. For the affine case,...
Definition: CoordinateElement.h:184
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition: CoordinateElement.h:109
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:42
void pull_back_nonaffine(mdspan2_t X, cmdspan2_t x, cmdspan2_t 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:62
void tabulate(int nd, std::span< const double > X, std::array< std::size_t, 2 > shape, std::span< double > basis) const
Evaluate basis values and derivatives at set of points.
Definition: CoordinateElement.cpp:47
bool is_affine() const noexcept
Check is geometry map is affine.
Definition: CoordinateElement.h:241
static double compute_jacobian_determinant(const U &J, std::span< typename U::value_type > w)
Compute the determinant of the Jacobian.
Definition: CoordinateElement.h:126
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell....
Definition: Expression.h:34
int value_size() const
Get value size.
Definition: Expression.h:231
void eval(const std::span< const std::int32_t > &cells, U &values) const
Evaluate the expression on cells.
Definition: Expression.h:140
std::pair< std::vector< double >, std::array< std::size_t, 2 > > X() const
Evaluation points on the reference cell.
Definition: Expression.h:243
std::shared_ptr< const FunctionSpace > argument_function_space() const
Get argument function space.
Definition: Expression.h:98
This class represents a function in a finite element function space , given by.
Definition: Function.h:45
void interpolate(const std::function< xt::xarray< T >(const xt::xtensor< double, 2 > &)> &f, std::span< const std::int32_t > cells)
Interpolate an expression function on a list of cells.
Definition: Function.h:174
void interpolate(const Expression< T > &e)
Interpolate an Expression (based on UFL) on all cells.
Definition: Function.h:301
void interpolate(const Function< T > &v)
Interpolate a Function.
Definition: Function.h:157
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(std::span< const double > x, std::array< std::size_t, 2 > xshape, std::span< const std::int32_t > cells, std::span< T > u, std::array< std::size_t, 2 > ushape) const
Evaluate the Function at points.
Definition: Function.h:326
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:222
std::shared_ptr< const FunctionSpace > function_space() const
Access the function space.
Definition: Function.h:136
std::string name
Name.
Definition: Function.h:577
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
void interpolate(const Function< T > &v, std::span< const std::int32_t > cells)
Interpolate a Function.
Definition: Function.h:150
void interpolate(const Expression< T > &e, std::span< const std::int32_t > cells)
Interpolate an Expression (based on UFL)
Definition: Function.h:242
~Function()=default
Destructor.
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:26
std::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:111
Distributed vector.
Definition: Vector.h:26
Finite element method functionality.
Definition: assemble_matrix_impl.h:25
void interpolate(Function< T > &u, std::span< const T > f, std::array< std::size_t, 2 > fshape, std::span< const std::int32_t > cells)
Interpolate an expression f(x) in a finite element space.
Definition: interpolate.h:421
std::vector< double > interpolation_coords(const FiniteElement &element, const mesh::Mesh &mesh, std::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:16