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>
27 void interpolate(Function<T>& u,
const Function<T>& v);
36 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
37 const Eigen::Ref<
const Eigen::Array<
double, 3, Eigen::Dynamic,
38 Eigen::RowMajor>>&)>& f);
51 const std::function<
void(
53 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
54 const Eigen::Ref<
const Eigen::Array<
double, Eigen::Dynamic, 3,
55 Eigen::RowMajor>>&)>& f);
63 void interpolate_values(
65 const Eigen::Ref<
const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic,
66 Eigen::RowMajor>>& values)
68 Eigen::Matrix<T, Eigen::Dynamic, 1>& coefficients = u.x()->array();
69 coefficients = Eigen::Map<const Eigen::Array<T, Eigen::Dynamic, 1>>(
70 values.data(), coefficients.rows());
74 void interpolate_from_any(Function<T>& u,
const Function<T>& v)
76 assert(v.function_space());
77 const auto element = u.function_space()->element();
79 if (!v.function_space()->has_element(*element))
81 throw std::runtime_error(
"Restricting finite elements function in "
82 "different elements not supported.");
85 const auto mesh = u.function_space()->mesh();
87 assert(v.function_space()->mesh());
88 if (mesh->id() != v.function_space()->mesh()->id())
90 throw std::runtime_error(
91 "Interpolation on different meshes not supported (yet).");
93 const int tdim = mesh->topology().dim();
96 assert(v.function_space());
97 std::shared_ptr<const fem::DofMap> dofmap_v = v.function_space()->dofmap();
99 auto map = mesh->topology().index_map(tdim);
102 Eigen::Matrix<T, Eigen::Dynamic, 1>& expansion_coefficients = u.x()->array();
105 const auto dofmap_u = u.function_space()->dofmap();
106 const Eigen::Matrix<T, Eigen::Dynamic, 1>& v_array = v.x()->array();
107 const int num_cells = map->size_local() + map->num_ghosts();
108 for (
int c = 0; c < num_cells; ++c)
110 auto dofs_v = dofmap_v->cell_dofs(c);
111 auto cell_dofs = dofmap_u->cell_dofs(c);
112 assert(dofs_v.size() == cell_dofs.size());
113 for (Eigen::Index i = 0; i < dofs_v.size(); ++i)
114 expansion_coefficients[cell_dofs[i]] = v_array[dofs_v[i]];
121 template <
typename T>
130 element->value_rank() != rank_v)
132 throw std::runtime_error(
"Cannot interpolate function into function space. "
134 + std::to_string(rank_v)
135 +
") does not match rank of function space ("
136 + std::to_string(element->value_rank()) +
")");
140 for (
int i = 0; i < element->value_rank(); ++i)
142 if (
int v_dim = v.
function_space()->element()->value_dimension(i);
143 element->value_dimension(i) != v_dim)
145 throw std::runtime_error(
146 "Cannot interpolate function into function space. "
148 + std::to_string(i) +
" of function (" + std::to_string(v_dim)
149 +
") does not match dimension " + std::to_string(i)
150 +
" of function space(" + std::to_string(element->value_dimension(i))
155 detail::interpolate_from_any(u, v);
158 template <
typename T>
162 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
163 const Eigen::Ref<
const Eigen::Array<
double, 3, Eigen::Dynamic,
164 Eigen::RowMajor>>&)>& f)
170 const Eigen::Array<double, 3, Eigen::Dynamic, Eigen::RowMajor> x
172 ->tabulate_scalar_subspace_dof_coordinates()
174 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values
179 std::vector<int> vshape(element->value_rank(), 1);
180 for (std::size_t i = 0; i < vshape.size(); ++i)
181 vshape[i] = element->value_dimension(i);
182 const int value_size = std::accumulate(std::begin(vshape), std::end(vshape),
183 1, std::multiplies<>());
187 const int tdim = mesh->topology().dim();
188 if (element->family() ==
"Mixed")
192 std::shared_ptr<const fem::DofMap> dofmap = u.
function_space()->dofmap();
193 auto map = mesh->topology().index_map(tdim);
195 const int num_cells = map->size_local() + map->num_ghosts();
197 Eigen::Array<T, 1, Eigen::Dynamic> mixed_values(values.cols());
198 int value_offset = 0;
199 for (
int i = 0; i < element->num_sub_elements(); ++i)
201 const std::vector<int> component = {i};
203 std::shared_ptr<const fem::FiniteElement> sub_element
204 = element->extract_sub_element(component);
205 const int element_block_size = sub_element->block_size();
206 for (std::size_t cell = 0; cell < num_cells; ++cell)
208 const auto cell_dofs = sub_dofmap.
cell_dofs(cell);
209 const std::size_t scalar_dofs = cell_dofs.rows() / element_block_size;
210 for (std::size_t dof = 0; dof < scalar_dofs; ++dof)
211 for (
int b = 0; b < element_block_size; ++b)
212 mixed_values(cell_dofs[element_block_size * dof + b])
213 = values(value_offset + b, cell_dofs[element_block_size * dof]);
215 value_offset += element_block_size;
217 detail::interpolate_values<T>(u, mixed_values);
224 if (values.cols() == 1 and values.rows() != 1)
226 if (values.rows() != x.cols())
228 throw std::runtime_error(
"Number of computed values is not equal to the "
229 "number of evaluation points. (1)");
231 detail::interpolate_values<T>(u, values);
235 if (values.rows() != value_size)
236 throw std::runtime_error(
"Values shape is incorrect. (2)");
237 if (values.cols() != x.cols())
239 throw std::runtime_error(
"Number of computed values is not equal to the "
240 "number of evaluation points. (2)");
243 detail::interpolate_values<T>(u, values.transpose());
247 template <
typename T>
250 const std::function<
void(
252 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>,
253 const Eigen::Ref<
const Eigen::Array<
double, Eigen::Dynamic, 3,
254 Eigen::RowMajor>>&)>& f)
257 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor> x
263 std::vector<int> vshape(element->value_rank(), 1);
264 for (std::size_t i = 0; i < vshape.size(); ++i)
265 vshape[i] = element->value_dimension(i);
266 const int value_size = std::accumulate(std::begin(vshape), std::end(vshape),
267 1, std::multiplies<>());
268 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> values(
269 x.rows(), value_size);
272 detail::interpolate_values<T>(u, values);
Degree-of-freedom map.
Definition: DofMap.h:66
Eigen::Array< std::int32_t, Eigen::Dynamic, 1 >::ConstSegmentReturnType cell_dofs(int cell) const
Local-to-global mapping of dofs on a cell.
Definition: DofMap.h:95
DofMap extract_sub_dofmap(const std::vector< int > &component) const
Extract subdofmap component.
Definition: DofMap.cpp:207
This class represents a function in a finite element function space , given by.
Definition: Function.h:42
std::shared_ptr< const FunctionSpace > function_space() const
Return shared pointer to function space.
Definition: Function.h:136
Functions tools, including FEM functions and pointwise defined functions.
Definition: assembler.h:19
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). This interface uses an expression function f that has an in/out argum...
Definition: interpolate.h:248
void interpolate(Function< T > &u, const Function< T > &v)
Interpolate a Function (on possibly non-matching meshes)
Definition: interpolate.h:122