9 #include "FunctionSpace.h"
10 #include <dolfinx/common/span.hpp>
11 #include <dolfinx/fem/DofMap.h>
12 #include <dolfinx/fem/FiniteElement.h>
13 #include <dolfinx/mesh/Mesh.h>
16 #include <xtensor/xtensor.hpp>
34 xt::xtensor<double, 2>
36 const tcb::span<const std::int32_t>& cells);
43 void interpolate(Function<T>& u,
const Function<T>& v);
59 const std::function<xt::xarray<T>(
const xt::xtensor<double, 2>&)>& f,
60 const xt::xtensor<double, 2>& x,
61 const tcb::span<const std::int32_t>& cells);
82 const std::function<
void(xt::xarray<T>&,
const xt::xtensor<double, 2>&)>& f,
83 const xt::xtensor<double, 2>& x,
84 const tcb::span<const std::int32_t>& cells);
90 void interpolate_from_any(Function<T>& u,
const Function<T>& v)
92 assert(v.function_space());
93 const auto element = u.function_space()->element();
95 if (v.function_space()->element()->hash() != element->hash())
97 throw std::runtime_error(
"Restricting finite elements function in "
98 "different elements not supported.");
101 const auto mesh = u.function_space()->mesh();
103 assert(v.function_space()->mesh());
104 if (mesh->id() != v.function_space()->mesh()->id())
106 throw std::runtime_error(
107 "Interpolation on different meshes not supported (yet).");
109 const int tdim = mesh->topology().dim();
112 assert(v.function_space());
113 std::shared_ptr<const fem::DofMap> dofmap_v = v.function_space()->dofmap();
115 auto map = mesh->topology().index_map(tdim);
118 std::vector<T>& coeffs = u.x()->mutable_array();
121 const auto dofmap_u = u.function_space()->dofmap();
122 const std::vector<T>& v_array = v.x()->array();
123 const int num_cells = map->size_local() + map->num_ghosts();
124 const int bs = dofmap_v->bs();
125 assert(bs == dofmap_u->bs());
126 for (
int c = 0; c < num_cells; ++c)
128 tcb::span<const std::int32_t> dofs_v = dofmap_v->cell_dofs(c);
129 tcb::span<const std::int32_t> cell_dofs = dofmap_u->cell_dofs(c);
130 assert(dofs_v.size() == cell_dofs.size());
131 for (std::size_t i = 0; i < dofs_v.size(); ++i)
133 for (
int k = 0; k < bs; ++k)
134 coeffs[bs * cell_dofs[i] + k] = v_array[bs * dofs_v[i] + k];
142 template <
typename T>
146 const std::shared_ptr<const FiniteElement> element
152 element->value_rank() != rank_v)
154 throw std::runtime_error(
"Cannot interpolate function into function space. "
156 + std::to_string(rank_v)
157 +
") does not match rank of function space ("
158 + std::to_string(element->value_rank()) +
")");
162 for (
int i = 0; i < element->value_rank(); ++i)
164 if (
int v_dim = v.
function_space()->element()->value_dimension(i);
165 element->value_dimension(i) != v_dim)
167 throw std::runtime_error(
168 "Cannot interpolate function into function space. "
170 + std::to_string(i) +
" of function (" + std::to_string(v_dim)
171 +
") does not match dimension " + std::to_string(i)
172 +
" of function space(" + std::to_string(element->value_dimension(i))
177 detail::interpolate_from_any(u, v);
180 template <
typename T>
183 const std::function<xt::xarray<T>(
const xt::xtensor<double, 2>&)>& f,
184 const xt::xtensor<double, 2>& x,
const tcb::span<const std::int32_t>& cells)
186 const std::shared_ptr<const FiniteElement> element
189 const int element_bs = element->block_size();
190 if (
int num_sub = element->num_sub_elements();
191 num_sub > 0 and num_sub != element_bs)
193 throw std::runtime_error(
"Cannot directly interpolate a mixed space. "
194 "Interpolate into subspaces.");
202 const int gdim = mesh->geometry().dim();
203 const int tdim = mesh->topology().dim();
209 throw std::runtime_error(
210 "Interpolation into this space is not yet supported.");
212 mesh->topology_mutable().create_entity_permutations();
213 const std::vector<std::uint32_t>& cell_info
214 = mesh->topology().get_cell_permutation_info();
221 xt::xarray<T> values = f(x);
223 if (values.dimension() == 1)
225 if (element->value_size() != 1)
226 throw std::runtime_error(
"Interpolation data has the wrong shape.");
227 values.reshape({element->value_size(), x.shape(1)});
230 if (values.shape(0) != element->value_size())
231 throw std::runtime_error(
"Interpolation data has the wrong shape.");
233 if (values.shape(1) != cells.size() * X.
shape[0])
234 throw std::runtime_error(
"Interpolation data has the wrong shape.");
239 const int dofmap_bs = dofmap->bs();
242 const int num_scalar_dofs = element->space_dimension() / element_bs;
243 const int value_size = element->value_size() / element_bs;
245 std::vector<T>& coeffs = u.
x()->mutable_array();
246 std::vector<T> _coeffs(num_scalar_dofs);
250 if (element->interpolation_ident())
252 for (std::int32_t c : cells)
254 tcb::span<const std::int32_t> dofs = dofmap->cell_dofs(c);
255 for (
int k = 0; k < element_bs; ++k)
257 for (
int i = 0; i < num_scalar_dofs; ++i)
258 _coeffs[i] = values(k, c * num_scalar_dofs + i);
259 element->apply_inverse_transpose_dof_transformation(_coeffs.data(),
261 for (
int i = 0; i < num_scalar_dofs; ++i)
263 const int dof = i * element_bs + k;
264 std::div_t pos = std::div(dof, dofmap_bs);
265 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
277 = mesh->geometry().dofmap();
279 const int num_dofs_g = x_dofmap.
num_links(0);
284 std::vector<double> J(X.
shape[0] * gdim * tdim);
285 std::vector<double> detJ(X.
shape[0]);
286 std::vector<double> K(X.
shape[0] * tdim * gdim);
296 xt::xtensor<double, 4> tabulated_data = cmap.
tabulate(1, X);
297 for (std::int32_t c : cells)
299 auto x_dofs = x_dofmap.
links(c);
300 for (
int i = 0; i < num_dofs_g; ++i)
301 for (
int j = 0; j < gdim; ++j)
302 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
308 tcb::span<const std::int32_t> dofs = dofmap->cell_dofs(c);
309 for (
int k = 0; k < element_bs; ++k)
312 for (
int m = 0; m < value_size; ++m)
314 std::copy_n(&values(k * value_size + m, c * X.
shape[0]), X.
shape[0],
315 _vals.
row(m).begin());
319 element->map_pull_back(reference_data.
data(), _vals.
data(), J.data(),
320 detJ.data(), K.data(), gdim, value_size, 1,
323 element->interpolate(reference_data, tcb::make_span(_coeffs));
324 element->apply_inverse_transpose_dof_transformation(_coeffs.data(),
327 assert(_coeffs.size() == num_scalar_dofs);
330 for (
int i = 0; i < num_scalar_dofs; ++i)
332 const int dof = i * element_bs + k;
333 std::div_t pos = std::div(dof, dofmap_bs);
334 coeffs[dofmap_bs * dofs[pos.quot] + pos.rem] = _coeffs[i];
341 template <
typename T>
344 const std::function<
void(xt::xarray<T>&,
const xt::xtensor<double, 2>&)>& f,
345 const xt::xtensor<double, 2>& x,
const tcb::span<const std::int32_t>& cells)
347 const std::shared_ptr<const FiniteElement> element
350 std::vector<int> vshape(element->value_rank(), 1);
351 for (std::size_t i = 0; i < vshape.size(); ++i)
352 vshape[i] = element->value_dimension(i);
353 const std::size_t value_size = std::accumulate(
354 std::begin(vshape), std::end(vshape), 1, std::multiplies<>());
356 auto fn = [value_size, &f](
const xt::xtensor<double, 2>& x) {
357 xt::xarray<T> values = xt::empty<T>({value_size, x.shape(1)});
362 interpolate<T>(u, fn, x, cells);
This class provides a dynamic 2-dimensional row-wise array data structure.
Definition: array2d.h:21
std::array< size_type, 2 > shape
The shape of the array.
Definition: array2d.h:157
constexpr tcb::span< value_type > row(size_type i)
Access a row in the array.
Definition: array2d.h:116
constexpr value_type * data() noexcept
Get pointer to the first element of the underlying storage.
Definition: array2d.h:133
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:31
xt::xtensor< double, 4 > tabulate(int n, const array2d< double > &X) const
Tabulate shape functions up to n-th order derivative at points X in the reference geometry Note: Dyna...
Definition: CoordinateElement.cpp:248
void compute_jacobian_data(const xt::xtensor< double, 4 > &tabulated_data, const array2d< double > &X, const array2d< double > &cell_geometry, std::vector< double > &J, tcb::span< double > detJ, std::vector< double > &K) const
Compute J, K and detJ for a cell with given geometry, and the basis functions and first order derivat...
Definition: CoordinateElement.cpp:255
This class represents a function in a finite element function space , given by.
Definition: Function.h:46
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:154
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:195
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:68
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:151
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:141
Finite element method functionality.
Definition: assemble_matrix_impl.h:23
xt::xtensor< double, 2 > interpolation_coords(const fem::FiniteElement &element, const mesh::Mesh &mesh, const tcb::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.cpp:17
void interpolate_c(Function< T > &u, const std::function< void(xt::xarray< T > &, const xt::xtensor< double, 2 > &)> &f, const xt::xtensor< double, 2 > &x, const tcb::span< const std::int32_t > &cells)
Interpolate an expression f(x)
Definition: interpolate.h:342
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:143