90 auto mesh = V0.mesh();
94 throw std::runtime_error(
"Meshes must be the same.");
96 if (
mesh->geometry().dim() != 3)
97 throw std::runtime_error(
"Geometric must be equal to 3..");
98 if (
mesh->geometry().dim() !=
mesh->topology()->dim())
100 throw std::runtime_error(
101 "Geometric and topological dimensions must be equal.");
103 constexpr int gdim = 3;
106 std::shared_ptr<const FiniteElement<T>> e0 = V0.element();
108 if (e0->map_type() != basix::maps::type::covariantPiola)
110 throw std::runtime_error(
111 "Finite element for parent space must be covariant Piola.");
114 std::shared_ptr<const FiniteElement<T>> e1 = V1.
element();
116 if (e1->map_type() != basix::maps::type::contravariantPiola)
118 throw std::runtime_error(
119 "Finite element for target space must be contracovariant Piola.");
123 std::span<const std::uint32_t> cell_info;
124 if (e1->needs_dof_transformations() or e0->needs_dof_transformations())
126 mesh->topology_mutable()->create_entity_permutations();
127 cell_info = std::span(
mesh->topology()->get_cell_permutation_info());
131 auto dofmap0 = V0.dofmap();
133 auto dofmap1 = V1.
dofmap();
137 auto apply_dof_transformation0
139 auto apply_inverse_dof_transform1 = e1->template dof_transformation_fn<U>(
143 const std::size_t space_dim0 = e0->space_dimension();
144 const std::size_t space_dim1 = e1->space_dimension();
145 if (e0->reference_value_size() != 3)
146 throw std::runtime_error(
"Value size for parent space should be 3.");
147 if (e1->reference_value_size() != 3)
148 throw std::runtime_error(
"Value size for target space should be 3.");
151 const auto [X, Xshape] = e1->interpolation_points();
155 auto x_dofmap =
mesh->geometry().dofmap();
156 const std::size_t num_dofs_g = cmap.
dim();
157 std::span<const T> x_g =
mesh->geometry().x();
158 std::array<std::size_t, 4> Phi_g_shape = cmap.
tabulate_shape(1, Xshape[0]);
159 std::vector<T> Phi_g_b(std::reduce(Phi_g_shape.begin(), Phi_g_shape.end(), 1,
162 md::mdspan<
const T, md::extents<std::size_t, gdim + 1, md::dynamic_extent,
163 md::dynamic_extent, 1>>
164 Phi_g(Phi_g_b.data(), Phi_g_shape);
165 cmap.
tabulate(1, X, Xshape, Phi_g_b);
168 std::vector<T> coord_dofs_b(num_dofs_g * gdim);
169 md::mdspan<T, md::extents<std::size_t, md::dynamic_extent, 3>> coord_dofs(
170 coord_dofs_b.data(), num_dofs_g, gdim);
171 std::vector<T> J_b(Xshape[0] * gdim * gdim);
172 md::mdspan<T, md::extents<std::size_t, md::dynamic_extent, 3, 3>> J(
173 J_b.data(), Xshape[0], gdim, gdim);
174 std::vector<T> K_b(Xshape[0] * gdim * gdim);
175 md::mdspan<T,
typename decltype(J)::extents_type> K(K_b.data(), J.extents());
176 std::vector<T> detJ(Xshape[0]);
177 std::vector<T> det_scratch(2 * gdim * gdim);
181 const auto [Phi0_b, Phi0_shape] = e0->tabulate(X, Xshape, 1);
182 md::mdspan<
const T, md::extents<std::size_t, 4, md::dynamic_extent,
183 md::dynamic_extent, 3>>
184 Phi0(Phi0_b.data(), Phi0_shape);
187 md::extents<std::size_t, md::dynamic_extent, md::dynamic_extent, 3, 3>
188 dPhi0_ext(Phi0.extent(1), Phi0.extent(2), Phi0.extent(0) - 1,
190 std::vector<T> dPhi0_b(dPhi0_ext.extent(0) * dPhi0_ext.extent(1)
191 * dPhi0_ext.extent(2) * dPhi0_ext.extent(3));
192 md::mdspan<T,
decltype(dPhi0_ext)> dPhi0(dPhi0_b.data(), dPhi0_ext);
197 const auto [Pi1_b, pi_shape] = e1->interpolation_operator();
198 md::mdspan<const T, md::dextents<std::size_t, 2>> Pi_1(Pi1_b.data(),
202 std::vector<T> curl_b(dPhi0.extent(0) * dPhi0.extent(1) * dPhi0.extent(3));
204 T, md::extents<std::size_t, md::dynamic_extent, md::dynamic_extent, 3>>
205 curl(curl_b.data(), dPhi0.extent(0), dPhi0.extent(1), dPhi0.extent(3));
207 std::vector<U> Ab(space_dim0 * space_dim1);
208 std::vector<U> local1(space_dim1);
211 assert(
mesh->topology()->index_map(gdim));
212 for (std::int32_t c = 0; c <
mesh->topology()->index_map(gdim)->size_local();
216 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
217 for (std::size_t i = 0; i < x_dofs.size(); ++i)
218 for (std::size_t j = 0; j < gdim; ++j)
219 coord_dofs(i, j) = x_g[3 * x_dofs[i] + j];
222 std::ranges::fill(J_b, 0);
223 for (std::size_t p = 0; p < Xshape[0]; ++p)
226 = md::submdspan(Phi_g, std::pair(1, gdim + 1), p, md::full_extent, 0);
227 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
229 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
240 for (std::size_t p = 0; p < Phi0.extent(1); ++p)
241 for (std::size_t phi = 0; phi < Phi0.extent(2); ++phi)
242 for (std::size_t d = 0; d < Phi0.extent(3); ++d)
243 for (std::size_t dx = 0; dx < dPhi0.extent(2); ++dx)
244 dPhi0(p, phi, dx, d) = Phi0(dx + 1, p, phi, d);
246 for (std::size_t p = 0; p < dPhi0.extent(0); ++p)
249 std::size_t size = dPhi0.extent(1) * dPhi0.extent(2) * dPhi0.extent(3);
250 std::size_t offset = p * size;
253 apply_dof_transformation0(std::span(dPhi0.data_handle() + offset, size),
255 dPhi0.extent(2) * dPhi0.extent(3));
261 for (std::size_t p = 0; p < curl.extent(0); ++p)
263 for (std::size_t i = 0; i < curl.extent(1); ++i)
265 curl(p, i, 0) = dPhi0(p, i, 1, 2) - dPhi0(p, i, 2, 1);
266 curl(p, i, 1) = dPhi0(p, i, 2, 0) - dPhi0(p, i, 0, 2);
267 curl(p, i, 2) = dPhi0(p, i, 0, 1) - dPhi0(p, i, 1, 0);
273 for (std::size_t i = 0; i < curl.extent(1); ++i)
275 auto values = md::submdspan(curl, md::full_extent, i, md::full_extent);
276 impl::interpolation_apply(Pi_1, values, std::span(local1), 1);
277 for (std::size_t j = 0; j < local1.size(); ++j)
278 Ab[space_dim0 * j + i] = local1[j];
281 apply_inverse_dof_transform1(Ab, cell_info, c, space_dim0);
282 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ab);
315 std::reference_wrapper<const DofMap>>
318 std::reference_wrapper<const DofMap>>
322 auto& e0 = V0.first.get();
323 const DofMap& dofmap0 = V0.second.get();
324 auto& e1 = V1.first.get();
325 const DofMap& dofmap1 = V1.second.get();
327 using cmdspan2_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
328 using cmdspan4_t = md::mdspan<const U, md::dextents<std::size_t, 4>>;
331 if (e0.map_type() != basix::maps::type::identity)
332 throw std::runtime_error(
"Wrong finite element space for V0.");
333 if (e0.block_size() != 1)
334 throw std::runtime_error(
"Block size is greater than 1 for V0.");
335 if (e0.reference_value_size() != 1)
336 throw std::runtime_error(
"Wrong value size for V0.");
338 if (e1.map_type() != basix::maps::type::covariantPiola)
339 throw std::runtime_error(
"Wrong finite element space for V1.");
340 if (e1.block_size() != 1)
341 throw std::runtime_error(
"Block size is greater than 1 for V1.");
344 const auto [X, Xshape] = e1.interpolation_points();
348 const int ndofs0 = e0.space_dimension();
349 const int tdim = topology.
dim();
350 std::vector<U> phi0_b((tdim + 1) * Xshape[0] * ndofs0 * 1);
351 cmdspan4_t phi0(phi0_b.data(), tdim + 1, Xshape[0], ndofs0, 1);
352 e0.tabulate(phi0_b, X, Xshape, 1);
356 cmdspan2_t dphi_reshaped(
357 phi0_b.data() + phi0.extent(3) * phi0.extent(2) * phi0.extent(1),
358 tdim * phi0.extent(1), phi0.extent(2));
361 auto apply_inverse_dof_transform = e1.template dof_transformation_fn<T>(
366 const std::vector<std::uint32_t>& cell_info
372 std::vector<T> Ab(e1.space_dimension() * ndofs0);
374 md::mdspan<T, md::dextents<std::size_t, 2>> A(Ab.data(),
375 e1.space_dimension(), ndofs0);
376 const auto [Pi, shape] = e1.interpolation_operator();
377 cmdspan2_t _Pi(Pi.data(), shape);
378 math::dot(_Pi, dphi_reshaped, A);
382 auto cell_map = topology.
index_map(tdim);
384 std::int32_t num_cells = cell_map->size_local();
385 std::vector<T> Ae(Ab.size());
386 for (std::int32_t c = 0; c < num_cells; ++c)
388 std::ranges::copy(Ab, Ae.begin());
389 apply_inverse_dof_transform(Ae, cell_info, c, ndofs0);
414 auto mesh = V0.mesh();
418 const int tdim =
mesh->topology()->dim();
419 const int gdim =
mesh->geometry().dim();
422 std::shared_ptr<const FiniteElement<U>> e0 = V0.element();
424 std::shared_ptr<const FiniteElement<U>> e1 = V1.
element();
427 std::span<const std::uint32_t> cell_info;
428 if (e1->needs_dof_transformations() or e0->needs_dof_transformations())
430 mesh->topology_mutable()->create_entity_permutations();
431 cell_info = std::span(
mesh->topology()->get_cell_permutation_info());
435 auto dofmap0 = V0.dofmap();
437 auto dofmap1 = V1.
dofmap();
441 const int bs0 = e0->block_size();
442 const int bs1 = e1->block_size();
443 auto apply_dof_transformation0
445 auto apply_inverse_dof_transform1 = e1->template dof_transformation_fn<T>(
449 const std::size_t space_dim0 = e0->space_dimension();
450 const std::size_t space_dim1 = e1->space_dimension();
451 const std::size_t dim0 = space_dim0 / bs0;
452 const std::size_t value_size_ref0 = e0->reference_value_size();
453 const std::size_t value_size0 = V0.element()->reference_value_size();
457 auto x_dofmap =
mesh->geometry().dofmap();
458 const std::size_t num_dofs_g = cmap.
dim();
459 std::span<const U> x_g =
mesh->geometry().x();
461 using mdspan2_t = md::mdspan<U, md::dextents<std::size_t, 2>>;
462 using cmdspan2_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
463 using cmdspan4_t = md::mdspan<const U, md::dextents<std::size_t, 4>>;
464 using mdspan3_t = md::mdspan<U, md::dextents<std::size_t, 3>>;
467 const auto [X, Xshape] = e1->interpolation_points();
468 std::array<std::size_t, 4> phi_shape = cmap.
tabulate_shape(1, Xshape[0]);
469 std::vector<U> phi_b(
470 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
471 cmdspan4_t phi(phi_b.data(), phi_shape);
475 std::vector<U> basis_derivatives_reference0_b(Xshape[0] * dim0
477 cmdspan4_t basis_derivatives_reference0(basis_derivatives_reference0_b.data(),
478 1, Xshape[0], dim0, value_size_ref0);
479 e0->tabulate(basis_derivatives_reference0_b, X, Xshape, 0);
482 std::ranges::transform(
483 basis_derivatives_reference0_b, basis_derivatives_reference0_b.begin(),
484 [atol = 1e-14](
auto x) { return std::abs(x) < atol ? 0.0 : x; });
487 std::vector<U> basis_reference0_b(Xshape[0] * dim0 * value_size_ref0);
488 mdspan3_t basis_reference0(basis_reference0_b.data(), Xshape[0], dim0,
490 std::vector<U> J_b(Xshape[0] * gdim * tdim);
491 mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
492 std::vector<U> K_b(Xshape[0] * tdim * gdim);
493 mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
494 std::vector<U> detJ(Xshape[0]);
495 std::vector<U> det_scratch(2 * tdim * gdim);
500 const auto [_Pi_1, pi_shape] = e1->interpolation_operator();
501 cmdspan2_t Pi_1(_Pi_1.data(), pi_shape);
503 bool interpolation_ident = e1->interpolation_ident();
505 using u_t = md::mdspan<U, md::dextents<std::size_t, 2>>;
506 using U_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
507 using J_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
508 using K_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
509 auto push_forward_fn0
510 = e0->basix_element().template map_fn<u_t, U_t, J_t, K_t>();
514 std::vector<U> basis_values_b(Xshape[0] * bs0 * dim0
516 mdspan3_t basis_values(basis_values_b.data(), Xshape[0], bs0 * dim0,
518 std::vector<U> mapped_values_b(Xshape[0] * bs0 * dim0
520 mdspan3_t mapped_values(mapped_values_b.data(), Xshape[0], bs0 * dim0,
524 = e1->basix_element().template map_fn<u_t, U_t, K_t, J_t>();
526 std::vector<U> coord_dofs_b(num_dofs_g * gdim);
527 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
528 std::vector<U> basis0_b(Xshape[0] * dim0 * value_size0);
529 mdspan3_t basis0(basis0_b.data(), Xshape[0], dim0, value_size0);
532 std::vector<T> Ab(space_dim0 * space_dim1);
533 std::vector<T> local1(space_dim1);
536 auto cell_map =
mesh->topology()->index_map(tdim);
538 std::int32_t num_cells = cell_map->size_local();
539 for (std::int32_t c = 0; c < num_cells; ++c)
542 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
543 for (std::size_t i = 0; i < x_dofs.size(); ++i)
545 for (
int j = 0; j < gdim; ++j)
546 coord_dofs(i, j) = x_g[3 * x_dofs[i] + j];
550 std::ranges::fill(J_b, 0);
551 for (std::size_t p = 0; p < Xshape[0]; ++p)
554 = md::submdspan(phi, std::pair(1, tdim + 1), p, md::full_extent, 0);
555 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
557 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
564 for (std::size_t k0 = 0; k0 < basis_reference0.extent(0); ++k0)
565 for (std::size_t k1 = 0; k1 < basis_reference0.extent(1); ++k1)
566 for (std::size_t k2 = 0; k2 < basis_reference0.extent(2); ++k2)
567 basis_reference0(k0, k1, k2)
568 = basis_derivatives_reference0(0, k0, k1, k2);
569 for (std::size_t p = 0; p < Xshape[0]; ++p)
571 apply_dof_transformation0(
572 std::span(basis_reference0.data_handle() + p * dim0 * value_size_ref0,
573 dim0 * value_size_ref0),
574 cell_info, c, value_size_ref0);
577 for (std::size_t p = 0; p < basis0.extent(0); ++p)
579 auto _u = md::submdspan(basis0, p, md::full_extent, md::full_extent);
580 auto _U = md::submdspan(basis_reference0, p, md::full_extent,
582 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
583 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
584 push_forward_fn0(_u, _U, _J, detJ[p], _K);
588 for (std::size_t p = 0; p < Xshape[0]; ++p)
589 for (std::size_t i = 0; i < dim0; ++i)
590 for (std::size_t j = 0; j < value_size0; ++j)
591 for (
int k = 0; k < bs0; ++k)
592 basis_values(p, i * bs0 + k, j * bs0 + k) = basis0(p, i, j);
595 for (std::size_t p = 0; p < basis_values.extent(0); ++p)
598 = md::submdspan(basis_values, p, md::full_extent, md::full_extent);
600 = md::submdspan(mapped_values, p, md::full_extent, md::full_extent);
601 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
602 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
603 pull_back_fn1(_U, _u, _K, 1.0 / detJ[p], _J);
608 if (interpolation_ident)
610 md::mdspan<T, md::dextents<std::size_t, 3>> A(
611 Ab.data(), Xshape[0], V1.
element()->value_size(), space_dim0);
612 for (std::size_t i = 0; i < mapped_values.extent(0); ++i)
613 for (std::size_t j = 0; j < mapped_values.extent(1); ++j)
614 for (std::size_t k = 0; k < mapped_values.extent(2); ++k)
615 A(i, k, j) = mapped_values(i, j, k);
619 for (std::size_t i = 0; i < mapped_values.extent(1); ++i)
622 = md::submdspan(mapped_values, md::full_extent, i, md::full_extent);
623 impl::interpolation_apply(Pi_1, values, std::span(local1), bs1);
624 for (std::size_t j = 0; j < local1.size(); j++)
625 Ab[space_dim0 * j + i] = local1[j];
629 apply_inverse_dof_transform1(Ab, cell_info, c, space_dim0);
630 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ab);