54 std::reference_wrapper<const DofMap>>
57 std::reference_wrapper<const DofMap>>
61 auto& e0 = V0.first.get();
62 const DofMap& dofmap0 = V0.second.get();
63 auto& e1 = V1.first.get();
64 const DofMap& dofmap1 = V1.second.get();
66 using cmdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
67 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
68 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
69 U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
70 using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
71 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
74 if (e0.map_type() != basix::maps::type::identity)
75 throw std::runtime_error(
"Wrong finite element space for V0.");
76 if (e0.block_size() != 1)
77 throw std::runtime_error(
"Block size is greater than 1 for V0.");
78 if (e0.reference_value_size() != 1)
79 throw std::runtime_error(
"Wrong value size for V0.");
81 if (e1.map_type() != basix::maps::type::covariantPiola)
82 throw std::runtime_error(
"Wrong finite element space for V1.");
83 if (e1.block_size() != 1)
84 throw std::runtime_error(
"Block size is greater than 1 for V1.");
87 const auto [X, Xshape] = e1.interpolation_points();
91 const int ndofs0 = e0.space_dimension();
92 const int tdim = topology.
dim();
93 std::vector<U> phi0_b((tdim + 1) * Xshape[0] * ndofs0 * 1);
94 cmdspan4_t phi0(phi0_b.data(), tdim + 1, Xshape[0], ndofs0, 1);
95 e0.tabulate(phi0_b, X, Xshape, 1);
99 cmdspan2_t dphi_reshaped(
100 phi0_b.data() + phi0.extent(3) * phi0.extent(2) * phi0.extent(1),
101 tdim * phi0.extent(1), phi0.extent(2));
104 auto apply_inverse_dof_transform
105 = e1.template get_pre_dof_transformation_function<T>(
110 const std::vector<std::uint32_t>& cell_info
116 std::vector<T> Ab(e1.space_dimension() * ndofs0);
118 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
119 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
120 A(Ab.data(), e1.space_dimension(), ndofs0);
121 const auto [Pi, shape] = e1.interpolation_operator();
122 cmdspan2_t _Pi(Pi.data(), shape);
123 math::dot(_Pi, dphi_reshaped, A);
127 auto cell_map = topology.
index_map(tdim);
129 std::int32_t num_cells = cell_map->size_local();
130 std::vector<T> Ae(Ab.size());
131 for (std::int32_t c = 0; c < num_cells; ++c)
133 std::copy(Ab.cbegin(), Ab.cend(), Ae.begin());
134 apply_inverse_dof_transform(Ae, cell_info, c, ndofs0);
159 auto mesh = V0.
mesh();
163 const int tdim = mesh->topology()->dim();
164 const int gdim = mesh->geometry().dim();
167 std::shared_ptr<const FiniteElement<U>> e0 = V0.
element();
169 std::shared_ptr<const FiniteElement<U>> e1 = V1.
element();
172 std::span<const std::uint32_t> cell_info;
173 if (e1->needs_dof_transformations() or e0->needs_dof_transformations())
175 mesh->topology_mutable()->create_entity_permutations();
176 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
180 auto dofmap0 = V0.
dofmap();
182 auto dofmap1 = V1.
dofmap();
186 const int bs0 = e0->block_size();
187 const int bs1 = e1->block_size();
188 auto apply_dof_transformation0
189 = e0->template get_pre_dof_transformation_function<U>(
191 auto apply_inverse_dof_transform1
192 = e1->template get_pre_dof_transformation_function<T>(
196 const std::size_t space_dim0 = e0->space_dimension();
197 const std::size_t space_dim1 = e1->space_dimension();
198 const std::size_t dim0 = space_dim0 / bs0;
199 const std::size_t value_size_ref0 = e0->reference_value_size() / bs0;
200 const std::size_t value_size0 = V0.
value_size() / bs0;
201 const std::size_t value_size1 = V1.
value_size() / bs1;
205 auto x_dofmap = mesh->geometry().dofmap();
206 const std::size_t num_dofs_g = cmap.
dim();
207 std::span<const U> x_g = mesh->geometry().x();
209 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
210 U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
211 using cmdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
212 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
213 using cmdspan3_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
214 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>;
215 using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
216 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
217 using mdspan3_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
218 U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>;
221 const auto [X, Xshape] = e1->interpolation_points();
222 std::array<std::size_t, 4> phi_shape = cmap.
tabulate_shape(1, Xshape[0]);
223 std::vector<U> phi_b(
224 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
225 cmdspan4_t phi(phi_b.data(), phi_shape);
229 std::vector<U> basis_derivatives_reference0_b(Xshape[0] * dim0
231 cmdspan4_t basis_derivatives_reference0(basis_derivatives_reference0_b.data(),
232 1, Xshape[0], dim0, value_size_ref0);
233 e0->tabulate(basis_derivatives_reference0_b, X, Xshape, 0);
236 std::transform(basis_derivatives_reference0_b.begin(),
237 basis_derivatives_reference0_b.end(),
238 basis_derivatives_reference0_b.begin(),
239 [atol = 1e-14](
auto x)
240 { return std::abs(x) < atol ? 0.0 : x; });
243 std::vector<U> basis_reference0_b(Xshape[0] * dim0 * value_size_ref0);
244 mdspan3_t basis_reference0(basis_reference0_b.data(), Xshape[0], dim0,
246 std::vector<U> J_b(Xshape[0] * gdim * tdim);
247 mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
248 std::vector<U> K_b(Xshape[0] * tdim * gdim);
249 mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
250 std::vector<U> detJ(Xshape[0]);
251 std::vector<U> det_scratch(2 * tdim * gdim);
256 const auto [_Pi_1, pi_shape] = e1->interpolation_operator();
257 cmdspan2_t Pi_1(_Pi_1.data(), pi_shape);
259 bool interpolation_ident = e1->interpolation_ident();
261 using u_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
262 U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
263 using U_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
264 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
265 using J_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
266 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
267 using K_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
268 const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
269 auto push_forward_fn0
270 = e0->basix_element().template map_fn<u_t, U_t, J_t, K_t>();
274 std::vector<U> basis_values_b(Xshape[0] * bs0 * dim0 * V1.
value_size());
275 mdspan3_t basis_values(basis_values_b.data(), Xshape[0], bs0 * dim0,
277 std::vector<U> mapped_values_b(Xshape[0] * bs0 * dim0 * V1.
value_size());
278 mdspan3_t mapped_values(mapped_values_b.data(), Xshape[0], bs0 * dim0,
282 = e1->basix_element().template map_fn<u_t, U_t, K_t, J_t>();
284 std::vector<U> coord_dofs_b(num_dofs_g * gdim);
285 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
286 std::vector<U> basis0_b(Xshape[0] * dim0 * value_size0);
287 mdspan3_t basis0(basis0_b.data(), Xshape[0], dim0, value_size0);
290 std::vector<T> Ab(space_dim0 * space_dim1);
291 std::vector<T> local1(space_dim1);
294 auto cell_map = mesh->topology()->index_map(tdim);
296 std::int32_t num_cells = cell_map->size_local();
297 for (std::int32_t c = 0; c < num_cells; ++c)
300 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
301 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
302 for (std::size_t i = 0; i < x_dofs.size(); ++i)
304 for (std::size_t j = 0; j < gdim; ++j)
305 coord_dofs(i, j) = x_g[3 * x_dofs[i] + j];
309 std::fill(J_b.begin(), J_b.end(), 0);
310 for (std::size_t p = 0; p < Xshape[0]; ++p)
312 auto dphi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
313 phi, std::pair(1, tdim + 1), p,
314 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
315 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
316 J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
317 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
319 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
320 K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
321 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
328 for (std::size_t k0 = 0; k0 < basis_reference0.extent(0); ++k0)
329 for (std::size_t k1 = 0; k1 < basis_reference0.extent(1); ++k1)
330 for (std::size_t k2 = 0; k2 < basis_reference0.extent(2); ++k2)
331 basis_reference0(k0, k1, k2)
332 = basis_derivatives_reference0(0, k0, k1, k2);
333 for (std::size_t p = 0; p < Xshape[0]; ++p)
335 apply_dof_transformation0(
336 std::span(basis_reference0.data_handle() + p * dim0 * value_size_ref0,
337 dim0 * value_size_ref0),
338 cell_info, c, value_size_ref0);
341 for (std::size_t p = 0; p < basis0.extent(0); ++p)
343 auto _u = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
344 basis0, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
345 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
346 auto _U = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
347 basis_reference0, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
348 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
349 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
350 K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
351 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
352 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
353 J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
354 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
355 push_forward_fn0(_u, _U, _J, detJ[p], _K);
359 for (std::size_t p = 0; p < Xshape[0]; ++p)
360 for (std::size_t i = 0; i < dim0; ++i)
361 for (std::size_t j = 0; j < value_size0; ++j)
362 for (
int k = 0; k < bs0; ++k)
363 basis_values(p, i * bs0 + k, j * bs0 + k) = basis0(p, i, j);
366 for (std::size_t p = 0; p < basis_values.extent(0); ++p)
368 auto _u = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
369 basis_values, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
370 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
371 auto _U = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
372 mapped_values, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
373 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
374 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
375 K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
376 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
377 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
378 J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
379 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
380 pull_back_fn1(_U, _u, _K, 1.0 / detJ[p], _J);
385 if (interpolation_ident)
387 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
388 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
389 A(Ab.data(), Xshape[0], V1.
value_size(), space_dim0);
390 for (std::size_t i = 0; i < mapped_values.extent(0); ++i)
391 for (std::size_t j = 0; j < mapped_values.extent(1); ++j)
392 for (std::size_t k = 0; k < mapped_values.extent(2); ++k)
393 A(i, k, j) = mapped_values(i, j, k);
397 for (std::size_t i = 0; i < mapped_values.extent(1); ++i)
399 auto values = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
400 mapped_values, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, i,
401 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
402 impl::interpolation_apply(Pi_1, values, std::span(local1), bs1);
403 for (std::size_t j = 0; j < local1.size(); j++)
404 Ab[space_dim0 * j + i] = local1[j];
408 apply_inverse_dof_transform1(Ab, cell_info, c, space_dim0);
409 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ab);