287 std::span<const std::int32_t> cells)
290 assert(_function_space);
291 assert(_function_space->element());
292 std::size_t value_size = e.value_size();
293 if (e.argument_function_space())
294 throw std::runtime_error(
"Cannot interpolate Expression with Argument");
296 if (value_size != _function_space->value_size())
298 throw std::runtime_error(
299 "Function value size not equal to Expression value size");
304 auto [X0, shape0] = e.X();
305 auto [X1, shape1] = _function_space->element()->interpolation_points();
306 if (shape0 != shape1)
308 throw std::runtime_error(
309 "Function element interpolation points has different shape to "
310 "Expression interpolation points");
313 for (std::size_t i = 0; i < X0.size(); ++i)
315 if (std::abs(X0[i] - X1[i]) > 1.0e-10)
317 throw std::runtime_error(
"Function element interpolation points not "
318 "equal to Expression interpolation points");
324 std::size_t num_cells = cells.size();
325 std::size_t num_points = e.X().second[0];
326 std::vector<value_type> fdata(num_cells * num_points * value_size);
327 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
329 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
330 f(fdata.data(), num_cells, num_points, value_size);
333 assert(_function_space->mesh());
334 e.eval(*_function_space->mesh(), cells, fdata,
335 {num_cells, num_points * value_size});
342 std::vector<value_type> fdata1(num_cells * num_points * value_size);
343 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
344 value_type, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
345 f1(fdata1.data(), value_size, num_cells, num_points);
346 for (std::size_t i = 0; i < f.extent(0); ++i)
347 for (std::size_t j = 0; j < f.extent(1); ++j)
348 for (std::size_t k = 0; k < f.extent(2); ++k)
349 f1(k, i, j) = f(i, j, k);
353 std::span<const value_type>(fdata1.data(), fdata1.size()),
354 {value_size, num_cells * num_points}, cells);
384 void eval(std::span<const geometry_type>
x, std::array<std::size_t, 2> xshape,
385 std::span<const std::int32_t> cells, std::span<value_type> u,
386 std::array<std::size_t, 2> ushape)
const
391 assert(
x.size() == xshape[0] * xshape[1]);
392 assert(u.size() == ushape[0] * ushape[1]);
397 if (xshape[0] != cells.size())
399 throw std::runtime_error(
400 "Number of points and number of cells must be equal.");
403 if (xshape[0] != ushape[0])
405 throw std::runtime_error(
406 "Length of array for Function values must be the "
407 "same as the number of points.");
411 assert(_function_space);
412 auto mesh = _function_space->mesh();
414 const std::size_t gdim = mesh->geometry().dim();
415 const std::size_t tdim = mesh->topology()->dim();
416 auto map = mesh->topology()->index_map(tdim);
422 auto x_dofmap = mesh->geometry().dofmap();
423 const std::size_t num_dofs_g = cmap.
dim();
424 auto x_g = mesh->geometry().x();
427 auto element = _function_space->element();
429 const int bs_element = element->block_size();
430 const std::size_t reference_value_size
431 = element->reference_value_size() / bs_element;
432 const std::size_t value_size = _function_space->value_size() / bs_element;
433 const std::size_t space_dimension = element->space_dimension() / bs_element;
437 const int num_sub_elements = element->num_sub_elements();
438 if (num_sub_elements > 1 and num_sub_elements != bs_element)
440 throw std::runtime_error(
"Function::eval is not supported for mixed "
441 "elements. Extract subspaces.");
445 std::vector<value_type> coefficients(space_dimension * bs_element);
448 std::shared_ptr<const DofMap> dofmap = _function_space->dofmap();
450 const int bs_dof = dofmap->bs();
452 std::span<const std::uint32_t> cell_info;
453 if (element->needs_dof_transformations())
455 mesh->topology_mutable()->create_entity_permutations();
456 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
459 std::vector<geometry_type> coord_dofs_b(num_dofs_g * gdim);
460 impl::mdspan_t<geometry_type, 2> coord_dofs(coord_dofs_b.data(), num_dofs_g,
462 std::vector<geometry_type> xp_b(1 * gdim);
463 impl::mdspan_t<geometry_type, 2> xp(xp_b.data(), 1, gdim);
466 std::fill(u.data(), u.data() + u.size(), 0.0);
467 std::span<const value_type> _v = _x->array();
471 std::array<std::size_t, 4> phi0_shape = cmap.
tabulate_shape(1, 1);
472 std::vector<geometry_type> phi0_b(std::reduce(
473 phi0_shape.begin(), phi0_shape.end(), 1, std::multiplies{}));
474 impl::mdspan_t<const geometry_type, 4> phi0(phi0_b.data(), phi0_shape);
475 cmap.
tabulate(1, std::vector<geometry_type>(tdim), {1, tdim}, phi0_b);
476 auto dphi0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
477 phi0, std::pair(1, tdim + 1), 0,
478 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
483 std::vector<geometry_type> phi_b(
484 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
485 impl::mdspan_t<const geometry_type, 4> phi(phi_b.data(), phi_shape);
486 auto dphi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
487 phi, std::pair(1, tdim + 1), 0,
488 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
491 std::vector<geometry_type> Xb(xshape[0] * tdim);
492 impl::mdspan_t<geometry_type, 2> X(Xb.data(), xshape[0], tdim);
495 std::vector<geometry_type> J_b(xshape[0] * gdim * tdim);
496 impl::mdspan_t<geometry_type, 3> J(J_b.data(), xshape[0], gdim, tdim);
497 std::vector<geometry_type> K_b(xshape[0] * tdim * gdim);
498 impl::mdspan_t<geometry_type, 3> K(K_b.data(), xshape[0], tdim, gdim);
499 std::vector<geometry_type> detJ(xshape[0]);
500 std::vector<geometry_type> det_scratch(2 * gdim * tdim);
503 for (std::size_t p = 0; p < cells.size(); ++p)
505 const int cell_index = cells[p];
512 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
513 x_dofmap, cell_index, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
514 assert(x_dofs.size() == num_dofs_g);
515 for (std::size_t i = 0; i < num_dofs_g; ++i)
517 const int pos = 3 * x_dofs[i];
518 for (std::size_t j = 0; j < gdim; ++j)
519 coord_dofs(i, j) = x_g[pos + j];
522 for (std::size_t j = 0; j < gdim; ++j)
523 xp(0, j) =
x[p * xshape[1] + j];
525 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
526 J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
527 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
528 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
529 K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
530 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
532 std::array<geometry_type, 3> Xpb = {0, 0, 0};
533 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
535 MDSPAN_IMPL_STANDARD_NAMESPACE::extents<
536 std::size_t, 1, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>
537 Xp(Xpb.data(), 1, tdim);
545 std::array<geometry_type, 3> x0 = {0, 0, 0};
546 for (std::size_t i = 0; i < coord_dofs.extent(1); ++i)
547 x0[i] += coord_dofs(0, i);
558 cmap.
tabulate(1, std::span(Xpb.data(), tdim), {1, tdim}, phi_b);
567 for (std::size_t j = 0; j < X.extent(1); ++j)
572 std::vector<geometry_type> basis_derivatives_reference_values_b(
573 1 * xshape[0] * space_dimension * reference_value_size);
574 impl::mdspan_t<const geometry_type, 4> basis_derivatives_reference_values(
575 basis_derivatives_reference_values_b.data(), 1, xshape[0],
576 space_dimension, reference_value_size);
577 std::vector<geometry_type> basis_values_b(space_dimension * value_size);
578 impl::mdspan_t<geometry_type, 2> basis_values(basis_values_b.data(),
579 space_dimension, value_size);
582 element->tabulate(basis_derivatives_reference_values_b, Xb,
583 {X.extent(0), X.extent(1)}, 0);
585 using xu_t = impl::mdspan_t<geometry_type, 2>;
586 using xU_t = impl::mdspan_t<const geometry_type, 2>;
587 using xJ_t = impl::mdspan_t<const geometry_type, 2>;
588 using xK_t = impl::mdspan_t<const geometry_type, 2>;
590 = element->basix_element().template map_fn<xu_t, xU_t, xJ_t, xK_t>();
592 auto apply_dof_transformation
594 ->template get_pre_dof_transformation_function<geometry_type>();
595 const std::size_t num_basis_values = space_dimension * reference_value_size;
597 for (std::size_t p = 0; p < cells.size(); ++p)
599 const int cell_index = cells[p];
607 apply_dof_transformation(
608 std::span(basis_derivatives_reference_values_b.data()
609 + p * num_basis_values,
611 cell_info, cell_index, reference_value_size);
614 auto _U = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
615 basis_derivatives_reference_values, 0, p,
616 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
617 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
618 auto _J = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
619 J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
620 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
621 auto _K = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
622 K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
623 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
624 push_forward_fn(basis_values, _U, _J, detJ[p], _K);
628 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell_index);
629 for (std::size_t i = 0; i < dofs.size(); ++i)
630 for (
int k = 0; k < bs_dof; ++k)
631 coefficients[bs_dof * i + k] = _v[bs_dof * dofs[i] + k];
634 for (
int k = 0; k < bs_element; ++k)
636 for (std::size_t i = 0; i < space_dimension; ++i)
638 for (std::size_t j = 0; j < value_size; ++j)
640 u[p * ushape[1] + (j * bs_element + k)]
641 += coefficients[bs_element * i + k] * basis_values(i, j);