10 #include "FunctionSpace.h"
11 #include <Eigen/Dense>
12 #include <dolfinx/fem/DofMap.h>
13 #include <dolfinx/fem/FiniteElement.h>
14 #include <dolfinx/mesh/Mesh.h>
28 void interpolate(Function<T>& u,
const Function<T>& v);
37 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
38 const Eigen::Ref<
const Eigen::Array<
double, 3, Eigen::Dynamic,
39 Eigen::RowMajor>>&)>& f);
55 const std::function<
void(
57 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
58 const Eigen::Ref<
const Eigen::Array<
double, Eigen::Dynamic, 3,
59 Eigen::RowMajor>>&)>& f);
65 void interpolate_from_any(Function<T>& u,
const Function<T>& v)
67 assert(v.function_space());
68 const auto element = u.function_space()->element();
70 if (v.function_space()->element()->hash() != element->hash())
72 throw std::runtime_error(
"Restricting finite elements function in "
73 "different elements not supported.");
76 const auto mesh = u.function_space()->mesh();
78 assert(v.function_space()->mesh());
79 if (mesh->id() != v.function_space()->mesh()->id())
81 throw std::runtime_error(
82 "Interpolation on different meshes not supported (yet).");
84 const int tdim = mesh->topology().dim();
87 assert(v.function_space());
88 std::shared_ptr<const fem::DofMap> dofmap_v = v.function_space()->dofmap();
90 auto map = mesh->topology().index_map(tdim);
93 std::vector<T>& coeffs = u.x()->mutable_array();
96 const auto dofmap_u = u.function_space()->dofmap();
97 const std::vector<T>& v_array = v.x()->array();
98 const int num_cells = map->size_local() + map->num_ghosts();
99 const int bs = dofmap_v->bs();
100 assert(bs == dofmap_u->bs());
101 for (
int c = 0; c < num_cells; ++c)
103 tcb::span<const std::int32_t> dofs_v = dofmap_v->cell_dofs(c);
104 tcb::span<const std::int32_t> cell_dofs = dofmap_u->cell_dofs(c);
105 assert(dofs_v.size() == cell_dofs.size());
106 for (std::size_t i = 0; i < dofs_v.size(); ++i)
108 for (
int k = 0; k < bs; ++k)
109 coeffs[bs * cell_dofs[i] + k] = v_array[bs * dofs_v[i] + k];
117 template <
typename T>
126 element->value_rank() != rank_v)
128 throw std::runtime_error(
"Cannot interpolate function into function space. "
130 + std::to_string(rank_v)
131 +
") does not match rank of function space ("
132 + std::to_string(element->value_rank()) +
")");
136 for (
int i = 0; i < element->value_rank(); ++i)
138 if (
int v_dim = v.
function_space()->element()->value_dimension(i);
139 element->value_dimension(i) != v_dim)
141 throw std::runtime_error(
142 "Cannot interpolate function into function space. "
144 + std::to_string(i) +
" of function (" + std::to_string(v_dim)
145 +
") does not match dimension " + std::to_string(i)
146 +
" of function space(" + std::to_string(element->value_dimension(i))
151 detail::interpolate_from_any(u, v);
154 template <
typename T>
158 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
159 const Eigen::Ref<
const Eigen::Array<
double, 3, Eigen::Dynamic,
160 Eigen::RowMajor>>&)>& f)
162 using EigenMatrixRowXd
163 = Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
167 const int element_bs = element->block_size();
169 if (
int num_sub = element->num_sub_elements();
170 num_sub > 0 and num_sub != element_bs)
172 throw std::runtime_error(
"Cannot directly interpolate a mixed space. "
173 "Interpolate into subspaces.");
180 const int tdim = mesh->topology().dim();
181 const int gdim = mesh->geometry().dim();
182 auto cell_map = mesh->topology().index_map(tdim);
184 const int num_cells = cell_map->size_local() + cell_map->num_ghosts();
188 = mesh->geometry().dofmap();
189 const int num_dofs_g = x_dofmap.
num_links(0);
190 const EigenMatrixRowXd& x_g = mesh->geometry().x();
194 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>& X
195 = element->interpolation_points();
197 mesh->topology_mutable().create_entity_permutations();
198 const std::vector<std::uint32_t>& cell_info
199 = mesh->topology().get_cell_permutation_info();
203 EigenMatrixRowXd x_cell(X.rows(), gdim);
204 std::vector<double> x;
205 EigenMatrixRowXd coordinate_dofs(num_dofs_g, gdim);
206 for (
int c = 0; c < num_cells; ++c)
209 auto x_dofs = x_dofmap.
links(c);
210 for (
int i = 0; i < num_dofs_g; ++i)
211 coordinate_dofs.row(i) = x_g.row(x_dofs[i]).head(gdim);
215 x.insert(x.end(), x_cell.data(), x_cell.data() + x_cell.size());
220 Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor> _x
221 = Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor>::Zero(
223 for (
int i = 0; i < gdim; ++i)
226 = Eigen::Map<Eigen::ArrayXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(
227 x.data() + i, x.size() / gdim,
228 Eigen::InnerStride<Eigen::Dynamic>(gdim));
236 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values;
241 if (element->value_size() == 1 and values.rows() > 1)
242 values = values.transpose().eval();
247 const int dofmap_bs = dofmap->bs();
253 const int num_scalar_dofs = element->space_dimension() / element_bs;
254 const int value_size = element->value_size() / element_bs;
257 if ((values.rows() != value_size * element_bs)
258 || (values.cols() != num_cells * X.rows()))
259 throw std::runtime_error(
"Interpolation data has the wrong shape.");
261 std::vector<T>& coeffs = u.
x()->mutable_array();
262 std::vector<T> _coeffs(num_scalar_dofs);
263 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> _vals;
264 for (
int c = 0; c < num_cells; ++c)
266 auto dofs = dofmap->cell_dofs(c);
267 for (
int k = 0; k < element_bs; ++k)
270 _vals = values.block(k * value_size, c * X.rows(), value_size, X.rows());
273 element->interpolate(_vals, cell_info[c], _coeffs);
274 assert(_coeffs.size() == num_scalar_dofs);
277 for (
int i = 0; i < num_scalar_dofs; ++i)
279 const int dof = i * element_bs + k;
280 std::div_t pos = std::div(dof, dofmap_bs);
281 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
287 template <
typename T>
290 const std::function<
void(
292 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
293 const Eigen::Ref<
const Eigen::Array<
double, Eigen::Dynamic, 3,
294 Eigen::RowMajor>>&)>& f)
298 std::vector<int> vshape(element->value_rank(), 1);
299 for (std::size_t i = 0; i < vshape.size(); ++i)
300 vshape[i] = element->value_dimension(i);
301 const int value_size = std::accumulate(std::begin(vshape), std::end(vshape),
302 1, std::multiplies<>());
306 &f](
const Eigen::Ref<
307 const Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor>>& x) {
308 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values(
309 x.cols(), value_size);
310 f(values, x.transpose());
314 interpolate<T>(u, fn);
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:26
void push_forward(Eigen::Ref< Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> x, 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 physical coordinates x for points X in the reference configuration.
Definition: CoordinateElement.cpp:67
This class represents a function in a finite element function space , given by.
Definition: Function.h:43
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:151
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:192
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
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
void interpolate_c(Function< T > &u, const std::function< void(Eigen::Ref< Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >>, const Eigen::Ref< const Eigen::Array< double, Eigen::Dynamic, 3, Eigen::RowMajor >> &)> &f)
Interpolate an expression f(x)
Definition: interpolate.h:288
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