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;
262 const int num_sub_elements = element->num_sub_elements();
263 if (num_sub_elements > 1 && num_sub_elements != block_size)
266 throw std::runtime_error(
267 "Blocked elements of mixed spaces are not yet supported.");
269 for (
int sub_e = 0; sub_e < num_sub_elements; ++sub_e)
271 std::shared_ptr<const fem::FiniteElement> sub_element
272 = element->extract_sub_element({sub_e});
274 const int sub_value_size = sub_element->value_size();
276 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> sub_u(
277 u.rows(), sub_value_size);
278 sub_f.eval(
x, cells, sub_u);
280 for (
int i = 0; i < sub_value_size; ++i)
281 u.col(offset + i) = sub_u.col(i);
282 offset += sub_value_size;
287 Eigen::Tensor<double, 3, Eigen::RowMajor> J(1, gdim, tdim);
288 Eigen::Array<double, Eigen::Dynamic, 1> detJ(1);
289 Eigen::Tensor<double, 3, Eigen::RowMajor> K(1, tdim, gdim);
290 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> X(
294 Eigen::Tensor<double, 3, Eigen::RowMajor> basis_reference_values(
295 1, space_dimension, reference_value_size);
296 Eigen::Tensor<double, 3, Eigen::RowMajor> basis_values(1, space_dimension,
300 Eigen::Matrix<T, 1, Eigen::Dynamic> coefficients(space_dimension
304 std::shared_ptr<const fem::DofMap> dofmap = _function_space->dofmap();
307 mesh->topology_mutable().create_entity_permutations();
308 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
309 = mesh->topology().get_cell_permutation_info();
313 const Eigen::Matrix<T, Eigen::Dynamic, 1>& _v = _x->array();
314 for (Eigen::Index p = 0; p < cells.rows(); ++p)
316 const int cell_index = cells(p);
323 auto x_dofs = x_dofmap.
links(cell_index);
324 for (
int i = 0; i < num_dofs_g; ++i)
325 coordinate_dofs.row(i) = x_g.row(x_dofs[i]).head(gdim);
332 element->evaluate_reference_basis(basis_reference_values, X);
335 element->transform_reference_basis(basis_values, basis_reference_values,
336 X, J, detJ, K, cell_info[cell_index]);
339 auto dofs = dofmap->cell_dofs(cell_index);
340 for (Eigen::Index i = 0; i < dofs.size(); ++i)
341 coefficients[i] = _v[dofs[i]];
344 for (
int block = 0; block < block_size; ++block)
346 for (
int i = 0; i < space_dimension; ++i)
348 for (
int j = 0; j < value_size; ++j)
351 u.row(p)[j * block_size + block]
352 += coefficients[i * block_size + block] * basis_values(0, i, j);
361 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
364 assert(_function_space);
365 std::shared_ptr<const mesh::Mesh> mesh = _function_space->mesh();
367 const int tdim = mesh->topology().dim();
370 const int value_size_loc = _function_space->element()->value_size();
373 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
374 point_values(mesh->geometry().x().rows(), value_size_loc);
378 = mesh->geometry().dofmap();
381 const int num_dofs_g = x_dofmap.
num_links(0);
382 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
383 = mesh->geometry().x();
387 Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>
x(num_dofs_g, 3);
388 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values(
389 num_dofs_g, value_size_loc);
390 auto map = mesh->topology().index_map(tdim);
392 const std::int32_t num_cells = map->size_local() + map->num_ghosts();
393 for (std::int32_t c = 0; c < num_cells; ++c)
396 auto dofs = x_dofmap.
links(c);
397 for (
int i = 0; i < num_dofs_g; ++i)
398 x.row(i) = x_g.row(dofs[i]);
400 values.resize(
x.rows(), value_size_loc);
403 Eigen::Array<int, Eigen::Dynamic, 1> cells(
x.rows());
405 eval(
x, cells, values);
408 for (
int i = 0; i <
x.rows(); ++i)
409 point_values.row(dofs[i]) = values.row(i);
419 std::size_t
id()
const {
return _id; }
426 std::shared_ptr<const FunctionSpace> _function_space;
429 std::shared_ptr<la::Vector<T>> _x;
432 mutable Vec _petsc_vector =
nullptr;