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>
47 explicit Function(std::shared_ptr<const FunctionSpace> V)
48 : _id(common::UniqueIdGenerator::
id()), _function_space(V),
49 _x(std::make_shared<la::Vector<T>>(V->dofmap()->index_map,
50 V->dofmap()->index_map_bs()))
52 if (!V->component().empty())
54 throw std::runtime_error(
"Cannot create Function from subspace. Consider "
55 "collapsing the function space");
65 Function(std::shared_ptr<const FunctionSpace> V,
67 : _id(common::UniqueIdGenerator::
id()), _function_space(V), _x(
x)
74 assert(V->dofmap()->index_map->size_global() * V->dofmap()->index_map_bs()
75 <= _x->bs() * _x->map()->size_global());
83 :
name(std::move(v.
name)), _id(std::move(v._id)),
84 _function_space(std::move(v._function_space)), _x(std::move(v._x)),
85 _petsc_vector(std::exchange(v._petsc_vector, nullptr))
93 VecDestroy(&_petsc_vector);
99 name = std::move(v.name);
100 _id = std::move(v._id);
101 _function_space = std::move(v._function_space);
102 _x = std::move(v._x);
103 std::swap(_petsc_vector, v._petsc_vector);
116 auto sub_space = _function_space->
sub({i});
127 const auto [function_space_new, collapsed_map]
131 assert(function_space_new);
132 auto vector_new = std::make_shared<la::Vector<T>>(
133 function_space_new->dofmap()->index_map,
134 function_space_new->dofmap()->index_map_bs());
137 const std::vector<T>& x_old = _x->array();
138 std::vector<T>& x_new = vector_new->mutable_array();
139 for (std::size_t i = 0; i < collapsed_map.size(); ++i)
141 assert((
int)i < x_new.size());
142 assert(collapsed_map[i] < x_old.size());
143 x_new[i] = x_old[collapsed_map[i]];
146 return Function(function_space_new, vector_new);
153 return _function_space;
162 assert(_function_space->dofmap());
163 assert(_function_space->dofmap()->index_map);
164 if (_x->bs() * _x->map()->size_global()
165 != _function_space->dofmap()->index_map->size_global()
166 * _function_space->dofmap()->index_map_bs())
168 throw std::runtime_error(
169 "Cannot access a non-const vector from a subfunction");
173 if constexpr (std::is_same<T, PetscScalar>::value)
178 *_function_space->dofmap()->index_map,
179 _function_space->dofmap()->index_map_bs(), _x->mutable_array());
182 return _petsc_vector;
186 throw std::runtime_error(
187 "Cannot return PETSc vector wrapper. Type mismatch");
192 std::shared_ptr<const la::Vector<T>>
x()
const {
return _x; }
195 std::shared_ptr<la::Vector<T>>
x() {
return _x; }
205 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
206 const Eigen::Ref<
const Eigen::Array<
double, 3, Eigen::Dynamic,
207 Eigen::RowMajor>>&)>& f)
224 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>>&
x,
225 const tcb::span<const std::int32_t>& cells,
227 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
233 if (
x.rows() != (
int)cells.size())
235 throw std::runtime_error(
236 "Number of points and number of cells must be equal.");
238 if (
x.rows() != u.rows())
240 throw std::runtime_error(
241 "Length of array for Function values must be the "
242 "same as the number of points.");
246 assert(_function_space);
247 std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
249 const int gdim = mesh->geometry().dim();
250 const int tdim = mesh->topology().dim();
251 auto map = mesh->topology().index_map(tdim);
255 = mesh->geometry().dofmap();
257 const int num_dofs_g = x_dofmap.
num_links(0);
258 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
259 = mesh->geometry().x();
265 assert(_function_space->element());
266 std::shared_ptr<const fem::FiniteElement> element
267 = _function_space->element();
269 const int bs_element = element->block_size();
270 const int reference_value_size
271 = element->reference_value_size() / bs_element;
272 const int value_size = element->value_size() / bs_element;
273 const int space_dimension = element->space_dimension() / bs_element;
277 const int num_sub_elements = element->num_sub_elements();
278 if (num_sub_elements > 1 and num_sub_elements != bs_element)
280 throw std::runtime_error(
"Function::eval is not supported for mixed "
281 "elements. Extract subspaces.");
285 Eigen::Tensor<double, 3, Eigen::RowMajor> J(1, gdim, tdim);
286 std::array<double, 1> detJ;
287 Eigen::Tensor<double, 3, Eigen::RowMajor> K(1, tdim, gdim);
288 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> X(
292 Eigen::Tensor<double, 3, Eigen::RowMajor> basis_reference_values(
293 1, space_dimension, reference_value_size);
294 Eigen::Tensor<double, 3, Eigen::RowMajor> basis_values(1, space_dimension,
298 Eigen::Matrix<T, 1, Eigen::Dynamic> coefficients(space_dimension
302 std::shared_ptr<const fem::DofMap> dofmap = _function_space->dofmap();
304 const int bs_dof = dofmap->bs();
306 mesh->topology_mutable().create_entity_permutations();
307 const std::vector<std::uint32_t>& cell_info
308 = mesh->topology().get_cell_permutation_info();
309 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
310 coordinate_dofs(num_dofs_g, gdim);
314 const std::vector<T>& _v = _x->mutable_array();
315 for (std::size_t p = 0; p < cells.size(); ++p)
317 const int cell_index = cells[p];
324 auto x_dofs = x_dofmap.
links(cell_index);
325 for (
int i = 0; i < num_dofs_g; ++i)
326 coordinate_dofs.row(i) = x_g.row(x_dofs[i]).head(gdim);
333 element->evaluate_reference_basis(basis_reference_values, X);
336 element->apply_dof_transformation(basis_reference_values.data(),
337 cell_info[cell_index],
338 reference_value_size);
341 element->transform_reference_basis(basis_values, basis_reference_values,
345 tcb::span<const std::int32_t> dofs = dofmap->cell_dofs(cell_index);
346 for (std::size_t i = 0; i < dofs.size(); ++i)
347 for (
int k = 0; k < bs_dof; ++k)
348 coefficients[bs_dof * i + k] = _v[bs_dof * dofs[i] + k];
351 auto u_row = u.row(p);
352 for (
int k = 0; k < bs_element; ++k)
354 for (
int i = 0; i < space_dimension; ++i)
356 for (
int j = 0; j < value_size; ++j)
359 u_row[j * bs_element + k]
360 += coefficients[bs_element * i + k] * basis_values(0, i, j);
369 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
372 assert(_function_space);
373 std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
375 const int tdim = mesh->topology().dim();
378 const int value_size_loc = _function_space->element()->value_size();
381 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
382 point_values(mesh->geometry().x().rows(), value_size_loc);
386 = mesh->geometry().dofmap();
389 const int num_dofs_g = x_dofmap.
num_links(0);
390 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
391 = mesh->geometry().x();
395 auto map = mesh->topology().index_map(tdim);
397 const std::int32_t num_cells = map->size_local() + map->num_ghosts();
399 std::vector<std::int32_t> cells(x_g.rows());
400 for (std::int32_t c = 0; c < num_cells; ++c)
403 tcb::span<const std::int32_t> dofs = x_dofmap.
links(c);
404 for (
int i = 0; i < num_dofs_g; ++i)
408 eval(x_g, cells, point_values);
417 std::size_t
id()
const {
return _id; }
424 std::shared_ptr<const FunctionSpace> _function_space;
427 std::shared_ptr<la::Vector<T>> _x;
430 mutable Vec _petsc_vector =
nullptr;
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:26
void compute_reference_geometry(Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &X, Eigen::Tensor< double, 3, Eigen::RowMajor > &J, tcb::span< double > 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) const
Compute reference coordinates X, and J, detJ and K for physical coordinates x.
Definition: CoordinateElement.cpp:85
This class represents a function in a finite element function space , given by.
Definition: Function.h:43
Vec vector() const
Return vector of expansion coefficients as a PETSc Vec. Throws an exception a PETSc Vec cannot be cre...
Definition: Function.h:159
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:65
void interpolate(const Function< T > &v)
Interpolate a Function (on possibly non-matching meshes)
Definition: Function.h:199
Function & operator=(Function &&v) noexcept
Move assignment.
Definition: Function.h:97
Function(Function &&v)
Move constructor.
Definition: Function.h:82
void eval(const Eigen::Ref< const Eigen::Array< double, Eigen::Dynamic, 3, Eigen::RowMajor >> &x, const tcb::span< const std::int32_t > &cells, Eigen::Ref< Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> u) const
Evaluate the Function at points.
Definition: Function.h:223
std::shared_ptr< la::Vector< T > > x()
Underlying vector.
Definition: Function.h:195
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:151
virtual ~Function()
Destructor.
Definition: Function.h:90
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:204
std::size_t id() const
ID.
Definition: Function.h:417
std::string name
Name.
Definition: Function.h:414
Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > compute_point_values() const
Compute values at all mesh 'nodes'.
Definition: Function.h:370
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:192
Function collapse() const
Collapse a subfunction (view into the Function) to a stand-alone Function.
Definition: Function.h:124
Function(std::shared_ptr< const FunctionSpace > V)
Create function on given function space.
Definition: Function.h:47
Function sub(int i) const
Extract subfunction (view into the Function)
Definition: Function.h:114
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:46
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:128
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:118
Distributed vector.
Definition: Vector.h:19
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
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:118
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:65