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>
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))
50 if (!V->component().empty())
52 throw std::runtime_error(
"Cannot create Function from subspace. Consider "
53 "collapsing the function space");
55 _x->array().setZero();
64 Function(std::shared_ptr<const FunctionSpace> V,
66 : _id(common::UniqueIdGenerator::
id()), _function_space(V), _x(
x)
73 assert(V->dofmap()->index_map->size_global()
74 * V->dofmap()->index_map->block_size()
75 <= _x->map()->block_size() * _x->map()->size_global());
88 VecDestroy(&_petsc_vector);
102 auto sub_space = _function_space->
sub({i});
113 const auto [function_space_new, collapsed_map]
117 assert(function_space_new);
118 auto vector_new = std::make_shared<la::Vector<T>>(
119 function_space_new->dofmap()->index_map);
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)
126 assert((
int)i < x_new.size());
127 assert(collapsed_map[i] < x_old.size());
128 x_new[i] = x_old[collapsed_map[i]];
131 return Function(function_space_new, vector_new);
138 return _function_space;
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())
153 throw std::runtime_error(
154 "Cannot access a non-const vector from a subfunction");
158 if constexpr (std::is_same<T, PetscScalar>::value)
163 *_function_space->dofmap()->index_map, _x->array());
165 return _petsc_vector;
169 throw std::runtime_error(
170 "Cannot return PETSc vector wrapper. Type mismatch");
175 std::shared_ptr<const la::Vector<T>>
x()
const {
return _x; }
178 std::shared_ptr<la::Vector<T>>
x() {
return _x; }
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)
207 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>>&
x,
208 const Eigen::Ref<
const Eigen::Array<int, Eigen::Dynamic, 1>>& cells,
210 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
216 if (
x.rows() != cells.rows())
218 throw std::runtime_error(
219 "Number of points and number of cells must be equal.");
221 if (
x.rows() != u.rows())
223 throw std::runtime_error(
224 "Length of array for Function values must be the "
225 "same as the number of points.");
229 assert(_function_space);
230 std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
232 const int gdim = mesh->geometry().dim();
233 const int tdim = mesh->topology().dim();
237 = mesh->geometry().dofmap();
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);
250 assert(_function_space->element());
251 std::shared_ptr<const fem::FiniteElement> element
252 = _function_space->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;
261 Eigen::Tensor<double, 3, Eigen::RowMajor> J(1, gdim, tdim);
262 Eigen::Array<double, Eigen::Dynamic, 1> detJ(1);
263 Eigen::Tensor<double, 3, Eigen::RowMajor> K(1, tdim, gdim);
264 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> X(
268 Eigen::Tensor<double, 3, Eigen::RowMajor> basis_reference_values(
269 1, space_dimension, reference_value_size);
270 Eigen::Tensor<double, 3, Eigen::RowMajor> basis_values(1, space_dimension,
274 Eigen::Matrix<T, 1, Eigen::Dynamic> coefficients(space_dimension
278 std::shared_ptr<const fem::DofMap> dofmap = _function_space->dofmap();
281 mesh->topology_mutable().create_entity_permutations();
282 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
283 = mesh->topology().get_cell_permutation_info();
287 const Eigen::Matrix<T, Eigen::Dynamic, 1>& _v = _x->array();
288 for (Eigen::Index p = 0; p < cells.rows(); ++p)
290 const int cell_index = cells(p);
297 auto x_dofs = x_dofmap.
links(cell_index);
298 for (
int i = 0; i < num_dofs_g; ++i)
299 coordinate_dofs.row(i) = x_g.row(x_dofs[i]).head(gdim);
306 element->evaluate_reference_basis(basis_reference_values, X);
309 element->transform_reference_basis(basis_values, basis_reference_values,
310 X, J, detJ, K, cell_info[cell_index]);
313 auto dofs = dofmap->cell_dofs(cell_index);
314 for (Eigen::Index i = 0; i < dofs.size(); ++i)
315 coefficients[i] = _v[dofs[i]];
318 for (
int block = 0; block < block_size; ++block)
320 for (
int i = 0; i < space_dimension; ++i)
322 for (
int j = 0; j < value_size; ++j)
325 u.row(p)[j * block_size + block]
326 += coefficients[i * block_size + block] * basis_values(0, i, j);
335 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
338 assert(_function_space);
339 std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
341 const int tdim = mesh->topology().dim();
344 const int value_size_loc = _function_space->element()->value_size();
347 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
348 point_values(mesh->geometry().x().rows(), value_size_loc);
352 = mesh->geometry().dofmap();
355 const int num_dofs_g = x_dofmap.
num_links(0);
356 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
357 = mesh->geometry().x();
361 Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>
x(num_dofs_g, 3);
362 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values(
363 num_dofs_g, value_size_loc);
364 auto map = mesh->topology().index_map(tdim);
366 const std::int32_t num_cells = map->size_local() + map->num_ghosts();
367 for (std::int32_t c = 0; c < num_cells; ++c)
370 auto dofs = x_dofmap.
links(c);
371 for (
int i = 0; i < num_dofs_g; ++i)
372 x.row(i) = x_g.row(dofs[i]);
374 values.resize(
x.rows(), value_size_loc);
377 Eigen::Array<int, Eigen::Dynamic, 1> cells(
x.rows());
379 eval(
x, cells, values);
382 for (
int i = 0; i <
x.rows(); ++i)
383 point_values.row(dofs[i]) = values.row(i);
393 std::size_t
id()
const {
return _id; }
400 std::shared_ptr<const FunctionSpace> _function_space;
403 std::shared_ptr<la::Vector<T>> _x;
406 mutable Vec _petsc_vector =
nullptr;