DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
Function.h
1// Copyright (C) 2003-2024 Anders Logg, Garth N. Wells and Massimiliano Leoni
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 "assembler.h"
13#include "interpolate.h"
14#include <algorithm>
15#include <basix/mdspan.hpp>
16#include <concepts>
17#include <dolfinx/common/IndexMap.h>
18#include <dolfinx/common/types.h>
19#include <dolfinx/la/Vector.h>
20#include <dolfinx/mesh/Geometry.h>
21#include <dolfinx/mesh/Mesh.h>
22#include <dolfinx/mesh/Topology.h>
23#include <functional>
24#include <memory>
25#include <numeric>
26#include <span>
27#include <string>
28#include <utility>
29#include <vector>
30
31namespace dolfinx::fem
32{
33template <dolfinx::scalar T, std::floating_point U>
34class Expression;
35
45template <dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_t<T>>
47{
48public:
51 using value_type = T;
53 using geometry_type = U;
54
57 explicit Function(std::shared_ptr<const FunctionSpace<geometry_type>> V)
58 : _function_space(V),
59 _x(std::make_shared<la::Vector<value_type>>(
60 V->dofmap()->index_map, V->dofmap()->index_map_bs()))
61 {
62 if (!V->component().empty())
63 {
64 throw std::runtime_error("Cannot create Function from subspace. Consider "
65 "collapsing the function space");
66 }
67 }
68
77 Function(std::shared_ptr<const FunctionSpace<geometry_type>> V,
78 std::shared_ptr<la::Vector<value_type>> x)
79 : _function_space(V), _x(x)
80 {
81 // NOTE: We do not check for a subspace since this constructor is
82 // used for creating subfunctions
83
84 // Assertion uses '<=' to deal with sub-functions
85 assert(V->dofmap());
86 assert(V->dofmap()->index_map->size_global() * V->dofmap()->index_map_bs()
87 <= _x->bs() * _x->index_map()->size_global());
88 }
89
90 // Copy constructor
91 Function(const Function& v) = delete;
92
94 Function(Function&& v) = default;
95
97 ~Function() = default;
98
100 Function& operator=(Function&& v) = default;
101
102 // Assignment
103 Function& operator=(const Function& v) = delete;
104
108 Function sub(int i) const
109 {
110 auto sub_space = std::make_shared<FunctionSpace<geometry_type>>(
111 _function_space->sub({i}));
112 assert(sub_space);
113 return Function(sub_space, _x);
114 }
115
120 {
121 // Create new collapsed FunctionSpace
122 auto [V, map] = _function_space->collapse();
123
124 // Create new vector
125 auto x = std::make_shared<la::Vector<value_type>>(
126 V.dofmap()->index_map, V.dofmap()->index_map_bs());
127
128 // Copy values into new vector
129 std::span<const value_type> x_old = _x->array();
130 std::span<value_type> x_new = x->mutable_array();
131 for (std::size_t i = 0; i < map.size(); ++i)
132 {
133 assert((int)i < x_new.size());
134 assert(map[i] < x_old.size());
135 x_new[i] = x_old[map[i]];
136 }
137
138 return Function(
139 std::make_shared<FunctionSpace<geometry_type>>(std::move(V)), x);
140 }
141
144 std::shared_ptr<const FunctionSpace<geometry_type>> function_space() const
145 {
146 return _function_space;
147 }
148
150 std::shared_ptr<const la::Vector<value_type>> x() const { return _x; }
151
153 std::shared_ptr<la::Vector<value_type>> x() { return _x; }
154
158 const std::function<
159 std::pair<std::vector<value_type>, std::vector<std::size_t>>(
160 md::mdspan<const geometry_type,
161 md::extents<std::size_t, 3, md::dynamic_extent>>)>& f)
162 {
163 assert(_function_space);
164 assert(_function_space->mesh());
165 int tdim = _function_space->mesh()->topology()->dim();
166 auto cmap = _function_space->mesh()->topology()->index_map(tdim);
167 assert(cmap);
168 std::vector<std::int32_t> cells(cmap->size_local() + cmap->num_ghosts(), 0);
169 std::iota(cells.begin(), cells.end(), 0);
170 interpolate(f, cells);
171 }
172
177 const std::function<
178 std::pair<std::vector<value_type>, std::vector<std::size_t>>(
179 md::mdspan<const geometry_type,
180 md::extents<std::size_t, 3, md::dynamic_extent>>)>& f,
181 std::span<const std::int32_t> cells)
182 {
183 assert(_function_space);
184 assert(_function_space->element());
185 assert(_function_space->mesh());
186 const std::vector<geometry_type> x
188 *_function_space->element(), _function_space->mesh()->geometry(),
189 cells);
190 md::mdspan<const geometry_type,
191 md::extents<std::size_t, 3, md::dynamic_extent>>
192 _x(x.data(), 3, x.size() / 3);
193
194 const auto [fx, fshape] = f(_x);
195 assert(fshape.size() <= 2);
196 if (int vs = _function_space->element()->value_size();
197 vs == 1 and fshape.size() == 1)
198 {
199 // Check for scalar-valued functions
200 if (fshape.front() != x.size() / 3)
201 throw std::runtime_error("Data returned by callable has wrong length");
202 }
203 else
204 {
205 // Check for vector/tensor value
206 if (fshape.size() != 2)
207 throw std::runtime_error("Expected 2D array of data");
208
209 if (fshape[0] != vs)
210 {
211 throw std::runtime_error(
212 "Data returned by callable has wrong shape(0) size");
213 }
214
215 if (fshape[1] != x.size() / 3)
216 {
217 throw std::runtime_error(
218 "Data returned by callable has wrong shape(1) size");
219 }
220 }
221
222 std::array<std::size_t, 2> _fshape;
223 if (fshape.size() == 1)
224 _fshape = {1, fshape[0]};
225 else
226 _fshape = {fshape[0], fshape[1]};
227
228 fem::interpolate(*this, std::span<const value_type>(fx.data(), fx.size()),
229 _fshape, cells);
230 }
231
238 {
239 assert(_function_space);
240 assert(_function_space->mesh());
241 int tdim = _function_space->mesh()->topology()->dim();
242 auto cmap = _function_space->mesh()->topology()->index_map(tdim);
243 assert(cmap);
244 std::vector<std::int32_t> cells(cmap->size_local() + cmap->num_ghosts(), 0);
245 std::iota(cells.begin(), cells.end(), 0);
246 interpolate(u, cells, cells);
247 }
248
266 std::span<const std::int32_t> cells0,
267 std::span<const std::int32_t> cells1 = {})
268 {
269 if (cells1.empty())
270 cells1 = cells0;
271 fem::interpolate(*this, cells1, u0, cells0);
272 }
273
280 {
281 assert(_function_space);
282 assert(_function_space->mesh());
283 int tdim = _function_space->mesh()->topology()->dim();
284 auto cmap = _function_space->mesh()->topology()->index_map(tdim);
285 assert(cmap);
286 std::vector<std::int32_t> cells(cmap->size_local() + cmap->num_ghosts(), 0);
287 std::iota(cells.begin(), cells.end(), 0);
288 interpolate(e, cells);
289 }
290
308 std::span<const std::int32_t> cells0,
309 std::span<const std::int32_t> cells1 = {})
310 {
311 // Extract mesh
312 const mesh::Mesh<geometry_type>* mesh0 = nullptr;
313 for (auto& c : e0.coefficients())
314 {
315 assert(c);
316 assert(c->function_space());
317 assert(c->function_space()->mesh());
319 = c->function_space()->mesh().get();
320 !mesh0)
321 {
322 mesh0 = mesh;
323 }
324 else if (mesh != mesh0)
325 {
326 throw std::runtime_error(
327 "Expression coefficient Functions have different meshes.");
328 }
329 }
330
331 // If Expression has no Function coefficients take mesh from `this`.
332 assert(_function_space);
333 assert(_function_space->mesh());
334 if (!mesh0)
335 mesh0 = _function_space->mesh().get();
336
337 // If cells1 is empty and Function and Expression meshes are the
338 // same, make cells1 the same as cells0. Otherwise check that
339 // lengths of cells0 and cells1 are the same.
340 if (cells1.empty() and mesh0 == _function_space->mesh().get())
341 cells1 = cells0;
342 else if (cells0.size() != cells1.size())
343 throw std::runtime_error("Cells lists have different lengths.");
344
345 // Check that Function and Expression spaces are compatible
346 assert(_function_space->element());
347 std::size_t value_size = e0.value_size();
348 if (e0.argument_space())
349 throw std::runtime_error("Cannot interpolate Expression with Argument.");
350 if (value_size != (std::size_t)_function_space->element()->value_size())
351 {
352 throw std::runtime_error(
353 "Function value size not equal to Expression value size.");
354 }
355
356 // Compatibility check
357 {
358 auto [X0, shape0] = e0.X();
359 auto [X1, shape1] = _function_space->element()->interpolation_points();
360 if (shape0 != shape1)
361 {
362 throw std::runtime_error(
363 "Function element interpolation points has different shape to "
364 "Expression interpolation points");
365 }
366
367 for (std::size_t i = 0; i < X0.size(); ++i)
368 {
369 if (std::abs(X0[i] - X1[i]) > 1.0e-10)
370 {
371 throw std::runtime_error("Function element interpolation points not "
372 "equal to Expression interpolation points");
373 }
374 }
375 }
376
377 // Array to hold evaluated Expression
378 std::size_t num_cells = cells0.size();
379 std::size_t num_points = e0.X().second[0];
380 std::vector<value_type> fdata(num_cells * num_points * value_size);
381 md::mdspan<const value_type, md::dextents<std::size_t, 3>> f(
382 fdata.data(), num_cells, num_points, value_size);
383
384 // Evaluate Expression at points
385 tabulate_expression(std::span(fdata), e0, *mesh0,
386 md::mdspan(cells0.data(), cells0.size()));
387
388 // Reshape evaluated data to fit interpolate.
389 // Expression returns matrix of shape (num_cells, num_points *
390 // value_size), i.e. xyzxyz ordering of dof values per cell per
391 // point. The interpolation uses xxyyzz input, ordered for all
392 // points of each cell, i.e. (value_size, num_cells*num_points).
393 std::vector<value_type> fdata1(num_cells * num_points * value_size);
394 md::mdspan<value_type, md::dextents<std::size_t, 3>> f1(
395 fdata1.data(), value_size, num_cells, num_points);
396 for (std::size_t i = 0; i < f.extent(0); ++i)
397 for (std::size_t j = 0; j < f.extent(1); ++j)
398 for (std::size_t k = 0; k < f.extent(2); ++k)
399 f1(k, i, j) = f(i, j, k);
400
401 // Interpolate values into appropriate space
402 fem::interpolate(*this,
403 std::span<const value_type>(fdata1.data(), fdata1.size()),
404 {value_size, num_cells * num_points}, cells1);
405 }
406
416 std::span<const std::int32_t> cells,
417 const geometry::PointOwnershipData<U>& interpolation_data)
418 {
419 fem::interpolate(*this, v, cells, interpolation_data);
420 }
421
434 void eval(std::span<const geometry_type> x, std::array<std::size_t, 2> xshape,
435 std::span<const std::int32_t> cells, std::span<value_type> u,
436 std::array<std::size_t, 2> ushape) const
437 {
438 if (cells.empty())
439 return;
440
441 assert(x.size() == xshape[0] * xshape[1]);
442 assert(u.size() == ushape[0] * ushape[1]);
443
444 // TODO: This could be easily made more efficient by exploiting
445 // points being ordered by the cell to which they belong.
446
447 if (xshape[0] != cells.size())
448 {
449 throw std::runtime_error(
450 "Number of points and number of cells must be equal.");
451 }
452
453 if (xshape[0] != ushape[0])
454 {
455 throw std::runtime_error(
456 "Length of array for Function values must be the "
457 "same as the number of points.");
458 }
459
460 // Get mesh
461 assert(_function_space);
462 auto mesh = _function_space->mesh();
463 assert(mesh);
464 const std::size_t gdim = mesh->geometry().dim();
465 const std::size_t tdim = mesh->topology()->dim();
466 auto map = mesh->topology()->index_map(tdim);
467
468 // Get coordinate map
469 const CoordinateElement<geometry_type>& cmap = mesh->geometry().cmap();
470
471 // Get geometry data
472 auto x_dofmap = mesh->geometry().dofmap();
473 const std::size_t num_dofs_g = cmap.dim();
474 auto x_g = mesh->geometry().x();
475
476 // Get element
477 auto element = _function_space->element();
478 assert(element);
479 const int bs_element = element->block_size();
480 const std::size_t reference_value_size = element->reference_value_size();
481 const std::size_t value_size
482 = _function_space->element()->reference_value_size();
483 const std::size_t space_dimension = element->space_dimension() / bs_element;
484
485 // If the space has sub elements, concatenate the evaluations on the
486 // sub elements
487 const int num_sub_elements = element->num_sub_elements();
488 if (num_sub_elements > 1 and num_sub_elements != bs_element)
489 {
490 throw std::runtime_error("Function::eval is not supported for mixed "
491 "elements. Extract subspaces.");
492 }
493
494 // Create work vector for expansion coefficients
495 std::vector<value_type> coefficients(space_dimension * bs_element);
496
497 // Get dofmap
498 std::shared_ptr<const DofMap> dofmap = _function_space->dofmap();
499 assert(dofmap);
500 const int bs_dof = dofmap->bs();
501
502 std::span<const std::uint32_t> cell_info;
503 if (element->needs_dof_transformations())
504 {
505 mesh->topology_mutable()->create_entity_permutations();
506 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
507 }
508
509 std::vector<geometry_type> coord_dofs_b(num_dofs_g * gdim);
510 impl::mdspan_t<geometry_type, 2> coord_dofs(coord_dofs_b.data(), num_dofs_g,
511 gdim);
512 std::vector<geometry_type> xp_b(1 * gdim);
513 impl::mdspan_t<geometry_type, 2> xp(xp_b.data(), 1, gdim);
514
515 // Loop over points
516 std::ranges::fill(u, 0.0);
517 std::span<const value_type> _v = _x->array();
518
519 // Evaluate geometry basis at point (0, 0, 0) on the reference cell.
520 // Used in affine case.
521 std::array<std::size_t, 4> phi0_shape = cmap.tabulate_shape(1, 1);
522 std::vector<geometry_type> phi0_b(std::reduce(
523 phi0_shape.begin(), phi0_shape.end(), 1, std::multiplies{}));
524 impl::mdspan_t<const geometry_type, 4> phi0(phi0_b.data(), phi0_shape);
525 cmap.tabulate(1, std::vector<geometry_type>(tdim), {1, tdim}, phi0_b);
526 auto dphi0
527 = md::submdspan(phi0, std::pair(1, tdim + 1), 0, md::full_extent, 0);
528
529 // Data structure for evaluating geometry basis at specific points.
530 // Used in non-affine case.
531 std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(1, 1);
532 std::vector<geometry_type> phi_b(
533 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
534 impl::mdspan_t<const geometry_type, 4> phi(phi_b.data(), phi_shape);
535 auto dphi
536 = md::submdspan(phi, std::pair(1, tdim + 1), 0, md::full_extent, 0);
537
538 // Reference coordinates for each point
539 std::vector<geometry_type> Xb(xshape[0] * tdim);
540 impl::mdspan_t<geometry_type, 2> X(Xb.data(), xshape[0], tdim);
541
542 // Geometry data at each point
543 std::vector<geometry_type> J_b(xshape[0] * gdim * tdim);
544 impl::mdspan_t<geometry_type, 3> J(J_b.data(), xshape[0], gdim, tdim);
545 std::vector<geometry_type> K_b(xshape[0] * tdim * gdim);
546 impl::mdspan_t<geometry_type, 3> K(K_b.data(), xshape[0], tdim, gdim);
547 std::vector<geometry_type> detJ(xshape[0]);
548 std::vector<geometry_type> det_scratch(2 * gdim * tdim);
549
550 // Prepare geometry data in each cell
551 for (std::size_t p = 0; p < cells.size(); ++p)
552 {
553 const int cell_index = cells[p];
554
555 // Skip negative cell indices
556 if (cell_index < 0)
557 continue;
558
559 // Get cell geometry (coordinate dofs)
560 auto x_dofs = md::submdspan(x_dofmap, cell_index, md::full_extent);
561 assert(x_dofs.size() == num_dofs_g);
562 for (std::size_t i = 0; i < num_dofs_g; ++i)
563 {
564 const int pos = 3 * x_dofs[i];
565 for (std::size_t j = 0; j < gdim; ++j)
566 coord_dofs(i, j) = x_g[pos + j];
567 }
568
569 for (std::size_t j = 0; j < gdim; ++j)
570 xp(0, j) = x[p * xshape[1] + j];
571
572 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
573 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
574
575 std::array<geometry_type, 3> Xpb = {0, 0, 0};
576 md::mdspan<geometry_type, md::extents<std::size_t, 1, md::dynamic_extent>>
577 Xp(Xpb.data(), 1, tdim);
578
579 // Compute reference coordinates X, and J, detJ and K
580 if (cmap.is_affine())
581 {
583 _J);
585 std::array<geometry_type, 3> x0 = {0, 0, 0};
586 for (std::size_t i = 0; i < coord_dofs.extent(1); ++i)
587 x0[i] += coord_dofs(0, i);
589 detJ[p]
591 _J, det_scratch);
592 }
593 else
594 {
595 // Pull-back physical point xp to reference coordinate Xp
596 cmap.pull_back_nonaffine(Xp, xp, coord_dofs);
597 cmap.tabulate(1, std::span(Xpb.data(), tdim), {1, tdim}, phi_b);
599 _J);
601 detJ[p]
603 _J, det_scratch);
604 }
605
606 for (std::size_t j = 0; j < X.extent(1); ++j)
607 X(p, j) = Xpb[j];
608 }
609
610 // Prepare basis function data structures
611 std::vector<geometry_type> basis_derivatives_reference_values_b(
612 1 * xshape[0] * space_dimension * reference_value_size);
613 impl::mdspan_t<const geometry_type, 4> basis_derivatives_reference_values(
614 basis_derivatives_reference_values_b.data(), 1, xshape[0],
615 space_dimension, reference_value_size);
616 std::vector<geometry_type> basis_values_b(space_dimension * value_size);
617 impl::mdspan_t<geometry_type, 2> basis_values(basis_values_b.data(),
618 space_dimension, value_size);
619
620 // Compute basis on reference element
621 element->tabulate(basis_derivatives_reference_values_b, Xb,
622 {X.extent(0), X.extent(1)}, 0);
623
624 using xu_t = impl::mdspan_t<geometry_type, 2>;
625 using xU_t = impl::mdspan_t<const geometry_type, 2>;
626 using xJ_t = impl::mdspan_t<const geometry_type, 2>;
627 using xK_t = impl::mdspan_t<const geometry_type, 2>;
628 auto push_forward_fn
629 = element->basix_element().template map_fn<xu_t, xU_t, xJ_t, xK_t>();
630
631 // Transformation function for basis function values
632 auto apply_dof_transformation
633 = element->template dof_transformation_fn<geometry_type>(
635
636 // Size of tensor for symmetric elements, unused in non-symmetric case, but
637 // placed outside the loop for pre-computation.
638 int matrix_size;
639 if (element->symmetric())
640 {
641 matrix_size = 0;
642 while (matrix_size * matrix_size < (int)ushape[1])
643 ++matrix_size;
644 }
645
646 const std::size_t num_basis_values = space_dimension * reference_value_size;
647 for (std::size_t p = 0; p < cells.size(); ++p)
648 {
649 const int cell_index = cells[p];
650 if (cell_index < 0) // Skip negative cell indices
651 continue;
652
653 // Permute the reference basis function values to account for the
654 // cell's orientation
655 apply_dof_transformation(
656 std::span(basis_derivatives_reference_values_b.data()
657 + p * num_basis_values,
658 num_basis_values),
659 cell_info, cell_index, reference_value_size);
660
661 {
662 auto _U = md::submdspan(basis_derivatives_reference_values, 0, p,
663 md::full_extent, md::full_extent);
664 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
665 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
666 push_forward_fn(basis_values, _U, _J, detJ[p], _K);
667 }
668
669 // Get degrees of freedom for current cell
670 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell_index);
671 for (std::size_t i = 0; i < dofs.size(); ++i)
672 for (int k = 0; k < bs_dof; ++k)
673 coefficients[bs_dof * i + k] = _v[bs_dof * dofs[i] + k];
674
675 if (element->symmetric())
676 {
677 int row = 0;
678 int rowstart = 0;
679 // Compute expansion
680 for (int k = 0; k < bs_element; ++k)
681 {
682 if (k - rowstart > row)
683 {
684 row++;
685 rowstart = k;
686 }
687 for (std::size_t i = 0; i < space_dimension; ++i)
688 {
689 for (std::size_t j = 0; j < value_size; ++j)
690 {
691 u[p * ushape[1]
692 + (j * bs_element + row * matrix_size + k - rowstart)]
693 += coefficients[bs_element * i + k] * basis_values(i, j);
694 if (k - rowstart != row)
695 {
696 u[p * ushape[1]
697 + (j * bs_element + row + matrix_size * (k - rowstart))]
698 += coefficients[bs_element * i + k] * basis_values(i, j);
699 }
700 }
701 }
702 }
703 }
704 else
705 {
706 // Compute expansion
707 for (int k = 0; k < bs_element; ++k)
708 {
709 for (std::size_t i = 0; i < space_dimension; ++i)
710 {
711 for (std::size_t j = 0; j < value_size; ++j)
712 {
713 u[p * ushape[1] + (j * bs_element + k)]
714 += coefficients[bs_element * i + k] * basis_values(i, j);
715 }
716 }
717 }
718 }
719 }
720 }
721
723 std::string name = "u";
724
725private:
726 // Function space
727 std::shared_ptr<const FunctionSpace<geometry_type>> _function_space;
728
729 // Vector of expansion coefficients (local)
730 std::shared_ptr<la::Vector<value_type>> _x;
731};
732
733} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
Functions supporting assembly of finite element fem::Form and fem::Expression.
Definition CoordinateElement.h:38
static void compute_jacobian(const U &dphi, const V &cell_geometry, W &&J)
Definition CoordinateElement.h:133
static void pull_back_affine(U &&X, const V &K, const std::array< T, 3 > &x0, const W &x)
Compute reference coordinates X for physical coordinates x for an affine map. For the affine case,...
Definition CoordinateElement.h:215
void tabulate(int nd, std::span< const T > X, std::array< std::size_t, 2 > shape, std::span< T > basis) const
Evaluate basis values and derivatives at set of points.
Definition CoordinateElement.cpp:55
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition CoordinateElement.h:142
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 tabulate.
Definition CoordinateElement.cpp:48
int dim() const
The dimension of the coordinate element space.
Definition CoordinateElement.cpp:205
void pull_back_nonaffine(mdspan2_t< T > X, mdspan2_t< const T > x, mdspan2_t< const T > cell_geometry, double tol=1.0e-6, int maxit=15) const
Compute reference coordinates X for physical coordinates x for a non-affine map.
Definition CoordinateElement.cpp:83
bool is_affine() const noexcept
Definition CoordinateElement.h:266
static double compute_jacobian_determinant(const U &J, std::span< typename U::value_type > w)
Compute the determinant of the Jacobian.
Definition CoordinateElement.h:159
An Expression represents a mathematical expression evaluated at a pre-defined points on a reference c...
Definition Expression.h:41
std::pair< std::vector< geometry_type >, std::array< std::size_t, 2 > > X() const
Evaluation point coordinates on the reference cell.
Definition Expression.h:166
const std::vector< std::shared_ptr< const Function< scalar_type, geometry_type > > > & coefficients() const
Expression coefficients.
Definition Expression.h:115
std::shared_ptr< const FunctionSpace< geometry_type > > argument_space() const
Argument function space.
Definition Expression.h:106
int value_size() const
Value size of the Expression result.
Definition Expression.h:156
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
Definition Function.h:47
U geometry_type
Geometry type of the Mesh that the Function is defined on.
Definition Function.h:53
void interpolate(const Function< value_type, geometry_type > &u)
Interpolate a Function over all cells.
Definition Function.h:237
std::shared_ptr< const FunctionSpace< geometry_type > > function_space() const
Access the function space.
Definition Function.h:144
void interpolate(const Function< value_type, geometry_type > &v, std::span< const std::int32_t > cells, const geometry::PointOwnershipData< U > &interpolation_data)
Interpolate a Function defined on a different mesh.
Definition Function.h:415
void interpolate(const Expression< value_type, geometry_type > &e0, std::span< const std::int32_t > cells0, std::span< const std::int32_t > cells1={})
Interpolate an Expression over a subset of cells.
Definition Function.h:307
Function(std::shared_ptr< const FunctionSpace< geometry_type > > V)
Create function on given function space.
Definition Function.h:57
void interpolate(const Function< value_type, geometry_type > &u0, std::span< const std::int32_t > cells0, std::span< const std::int32_t > cells1={})
Interpolate a Function over a subset of cells.
Definition Function.h:265
Function collapse() const
Collapse a subfunction (view into a Function) to a stand-alone Function.
Definition Function.h:119
void interpolate(const std::function< std::pair< std::vector< value_type >, std::vector< std::size_t > >(md::mdspan< const geometry_type, md::extents< std::size_t, 3, md::dynamic_extent > >)> &f)
Interpolate an expression f(x) on the whole domain.
Definition Function.h:157
Function sub(int i) const
Extract a sub-function (a view into the Function).
Definition Function.h:108
void interpolate(const std::function< std::pair< std::vector< value_type >, std::vector< std::size_t > >(md::mdspan< const geometry_type, md::extents< std::size_t, 3, md::dynamic_extent > >)> &f, std::span< const std::int32_t > cells)
Interpolate an expression f(x) over a set of cells.
Definition Function.h:176
std::string name
Name.
Definition Function.h:723
std::shared_ptr< la::Vector< value_type > > x()
Underlying vector.
Definition Function.h:153
void interpolate(const Expression< value_type, geometry_type > &e)
Interpolate an Expression on all cells.
Definition Function.h:279
~Function()=default
Destructor.
void eval(std::span< const geometry_type > x, std::array< std::size_t, 2 > xshape, std::span< const std::int32_t > cells, std::span< value_type > u, std::array< std::size_t, 2 > ushape) const
Evaluate the Function at points.
Definition Function.h:434
Function(Function &&v)=default
Move constructor.
std::shared_ptr< const la::Vector< value_type > > x() const
Underlying vector (const version).
Definition Function.h:150
Function & operator=(Function &&v)=default
Move assignment.
Function(std::shared_ptr< const FunctionSpace< geometry_type > > V, std::shared_ptr< la::Vector< value_type > > x)
Create function on given function space with a given vector.
Definition Function.h:77
T value_type
Definition Function.h:51
Definition Vector.h:32
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Finite element method functionality.
Definition assemble_expression_impl.h:23
void tabulate_expression(std::span< T > values, const fem::Expression< T, U > &e, md::mdspan< const T, md::dextents< std::size_t, 2 > > coeffs, std::span< const T > constants, const mesh::Mesh< U > &mesh, fem::MDSpan2 auto entities, std::optional< std::pair< std::reference_wrapper< const FiniteElement< U > >, std::size_t > > element)
Evaluate an Expression on cells or facets.
Definition assembler.h:64
void interpolate(Function< T, U > &u, std::span< const T > f, std::array< std::size_t, 2 > fshape, std::span< const std::int32_t > cells)
Interpolate an evaluated expression f(x) in a finite element space.
Definition interpolate.h:694
@ standard
Standard.
Definition FiniteElement.h:27
std::vector< T > interpolation_coords(const fem::FiniteElement< T > &element, const mesh::Geometry< T > &geometry, 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.h:51
Linear algebra interface.
Definition sparsitybuild.h:15
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
Information on the ownership of points distributed across processes.
Definition utils.h:30