308 std::span<const std::int32_t> cells0,
309 std::span<const std::int32_t> cells1 = {})
316 assert(c->function_space());
317 assert(c->function_space()->mesh());
319 = c->function_space()->mesh().get();
324 else if (
mesh != mesh0)
326 throw std::runtime_error(
327 "Expression coefficient Functions have different meshes.");
332 assert(_function_space);
333 assert(_function_space->mesh());
335 mesh0 = _function_space->mesh().get();
340 if (cells1.empty() and mesh0 == _function_space->mesh().get())
342 else if (cells0.size() != cells1.size())
343 throw std::runtime_error(
"Cells lists have different lengths.");
346 assert(_function_space->element());
349 throw std::runtime_error(
"Cannot interpolate Expression with Argument.");
350 if (value_size != (std::size_t)_function_space->element()->value_size())
352 throw std::runtime_error(
353 "Function value size not equal to Expression value size.");
358 auto [X0, shape0] = e0.
X();
359 auto [X1, shape1] = _function_space->element()->interpolation_points();
360 if (shape0 != shape1)
362 throw std::runtime_error(
363 "Function element interpolation points has different shape to "
364 "Expression interpolation points");
367 for (std::size_t i = 0; i < X0.size(); ++i)
369 if (std::abs(X0[i] - X1[i]) > 1.0e-10)
371 throw std::runtime_error(
"Function element interpolation points not "
372 "equal to Expression interpolation points");
378 std::size_t num_cells = cells0.size();
379 std::size_t num_points = e0.
X().second[0];
380 std::vector<value_type> fdata(num_cells * num_points * value_size);
381 md::mdspan<const value_type, md::dextents<std::size_t, 3>> f(
382 fdata.data(), num_cells, num_points, value_size);
386 md::mdspan(cells0.data(), cells0.size()));
393 std::vector<value_type> fdata1(num_cells * num_points * value_size);
394 md::mdspan<value_type, md::dextents<std::size_t, 3>> f1(
395 fdata1.data(), value_size, num_cells, num_points);
396 for (std::size_t i = 0; i < f.extent(0); ++i)
397 for (std::size_t j = 0; j < f.extent(1); ++j)
398 for (std::size_t k = 0; k < f.extent(2); ++k)
399 f1(k, i, j) = f(i, j, k);
403 std::span<const value_type>(fdata1.data(), fdata1.size()),
404 {value_size, num_cells * num_points}, cells1);
434 void eval(std::span<const geometry_type>
x, std::array<std::size_t, 2> xshape,
435 std::span<const std::int32_t> cells, std::span<value_type> u,
436 std::array<std::size_t, 2> ushape)
const
441 assert(
x.size() == xshape[0] * xshape[1]);
442 assert(u.size() == ushape[0] * ushape[1]);
447 if (xshape[0] != cells.size())
449 throw std::runtime_error(
450 "Number of points and number of cells must be equal.");
453 if (xshape[0] != ushape[0])
455 throw std::runtime_error(
456 "Length of array for Function values must be the "
457 "same as the number of points.");
461 assert(_function_space);
462 auto mesh = _function_space->mesh();
464 const std::size_t gdim =
mesh->geometry().dim();
465 const std::size_t tdim =
mesh->topology()->dim();
466 auto map =
mesh->topology()->index_map(tdim);
472 auto x_dofmap =
mesh->geometry().dofmap();
473 const std::size_t num_dofs_g = cmap.
dim();
474 auto x_g =
mesh->geometry().x();
477 auto element = _function_space->element();
479 const int bs_element = element->block_size();
480 const std::size_t reference_value_size = element->reference_value_size();
481 const std::size_t value_size
482 = _function_space->element()->reference_value_size();
483 const std::size_t space_dimension = element->space_dimension() / bs_element;
487 const int num_sub_elements = element->num_sub_elements();
488 if (num_sub_elements > 1 and num_sub_elements != bs_element)
490 throw std::runtime_error(
"Function::eval is not supported for mixed "
491 "elements. Extract subspaces.");
495 std::vector<value_type> coefficients(space_dimension * bs_element);
498 std::shared_ptr<const DofMap> dofmap = _function_space->dofmap();
500 const int bs_dof = dofmap->bs();
502 std::span<const std::uint32_t> cell_info;
503 if (element->needs_dof_transformations())
505 mesh->topology_mutable()->create_entity_permutations();
506 cell_info = std::span(
mesh->topology()->get_cell_permutation_info());
509 std::vector<geometry_type> coord_dofs_b(num_dofs_g * gdim);
510 impl::mdspan_t<geometry_type, 2> coord_dofs(coord_dofs_b.data(), num_dofs_g,
512 std::vector<geometry_type> xp_b(1 * gdim);
513 impl::mdspan_t<geometry_type, 2> xp(xp_b.data(), 1, gdim);
516 std::ranges::fill(u, 0.0);
517 std::span<const value_type> _v = _x->array();
521 std::array<std::size_t, 4> phi0_shape = cmap.
tabulate_shape(1, 1);
522 std::vector<geometry_type> phi0_b(std::reduce(
523 phi0_shape.begin(), phi0_shape.end(), 1, std::multiplies{}));
524 impl::mdspan_t<const geometry_type, 4> phi0(phi0_b.data(), phi0_shape);
525 cmap.
tabulate(1, std::vector<geometry_type>(tdim), {1, tdim}, phi0_b);
527 = md::submdspan(phi0, std::pair(1, tdim + 1), 0, md::full_extent, 0);
532 std::vector<geometry_type> phi_b(
533 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
534 impl::mdspan_t<const geometry_type, 4> phi(phi_b.data(), phi_shape);
536 = md::submdspan(phi, std::pair(1, tdim + 1), 0, md::full_extent, 0);
539 std::vector<geometry_type> Xb(xshape[0] * tdim);
540 impl::mdspan_t<geometry_type, 2> X(Xb.data(), xshape[0], tdim);
543 std::vector<geometry_type> J_b(xshape[0] * gdim * tdim);
544 impl::mdspan_t<geometry_type, 3> J(J_b.data(), xshape[0], gdim, tdim);
545 std::vector<geometry_type> K_b(xshape[0] * tdim * gdim);
546 impl::mdspan_t<geometry_type, 3> K(K_b.data(), xshape[0], tdim, gdim);
547 std::vector<geometry_type> detJ(xshape[0]);
548 std::vector<geometry_type> det_scratch(2 * gdim * tdim);
551 for (std::size_t p = 0; p < cells.size(); ++p)
553 const int cell_index = cells[p];
560 auto x_dofs = md::submdspan(x_dofmap, cell_index, md::full_extent);
561 assert(x_dofs.size() == num_dofs_g);
562 for (std::size_t i = 0; i < num_dofs_g; ++i)
564 const int pos = 3 * x_dofs[i];
565 for (std::size_t j = 0; j < gdim; ++j)
566 coord_dofs(i, j) = x_g[pos + j];
569 for (std::size_t j = 0; j < gdim; ++j)
570 xp(0, j) =
x[p * xshape[1] + j];
572 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
573 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
575 std::array<geometry_type, 3> Xpb = {0, 0, 0};
576 md::mdspan<geometry_type, md::extents<std::size_t, 1, md::dynamic_extent>>
577 Xp(Xpb.data(), 1, tdim);
585 std::array<geometry_type, 3> x0 = {0, 0, 0};
586 for (std::size_t i = 0; i < coord_dofs.extent(1); ++i)
587 x0[i] += coord_dofs(0, i);
597 cmap.
tabulate(1, std::span(Xpb.data(), tdim), {1, tdim}, phi_b);
606 for (std::size_t j = 0; j < X.extent(1); ++j)
611 std::vector<geometry_type> basis_derivatives_reference_values_b(
612 1 * xshape[0] * space_dimension * reference_value_size);
613 impl::mdspan_t<const geometry_type, 4> basis_derivatives_reference_values(
614 basis_derivatives_reference_values_b.data(), 1, xshape[0],
615 space_dimension, reference_value_size);
616 std::vector<geometry_type> basis_values_b(space_dimension * value_size);
617 impl::mdspan_t<geometry_type, 2> basis_values(basis_values_b.data(),
618 space_dimension, value_size);
621 element->tabulate(basis_derivatives_reference_values_b, Xb,
622 {X.extent(0), X.extent(1)}, 0);
624 using xu_t = impl::mdspan_t<geometry_type, 2>;
625 using xU_t = impl::mdspan_t<const geometry_type, 2>;
626 using xJ_t = impl::mdspan_t<const geometry_type, 2>;
627 using xK_t = impl::mdspan_t<const geometry_type, 2>;
629 = element->basix_element().template map_fn<xu_t, xU_t, xJ_t, xK_t>();
632 auto apply_dof_transformation
633 = element->template dof_transformation_fn<geometry_type>(
639 if (element->symmetric())
642 while (matrix_size * matrix_size < (
int)ushape[1])
646 const std::size_t num_basis_values = space_dimension * reference_value_size;
647 for (std::size_t p = 0; p < cells.size(); ++p)
649 const int cell_index = cells[p];
655 apply_dof_transformation(
656 std::span(basis_derivatives_reference_values_b.data()
657 + p * num_basis_values,
659 cell_info, cell_index, reference_value_size);
662 auto _U = md::submdspan(basis_derivatives_reference_values, 0, p,
663 md::full_extent, md::full_extent);
664 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
665 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
666 push_forward_fn(basis_values, _U, _J, detJ[p], _K);
670 std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell_index);
671 for (std::size_t i = 0; i < dofs.size(); ++i)
672 for (
int k = 0; k < bs_dof; ++k)
673 coefficients[bs_dof * i + k] = _v[bs_dof * dofs[i] + k];
675 if (element->symmetric())
680 for (
int k = 0; k < bs_element; ++k)
682 if (k - rowstart > row)
687 for (std::size_t i = 0; i < space_dimension; ++i)
689 for (std::size_t j = 0; j < value_size; ++j)
692 + (j * bs_element + row * matrix_size + k - rowstart)]
693 += coefficients[bs_element * i + k] * basis_values(i, j);
694 if (k - rowstart != row)
697 + (j * bs_element + row + matrix_size * (k - rowstart))]
698 += coefficients[bs_element * i + k] * basis_values(i, j);
707 for (
int k = 0; k < bs_element; ++k)
709 for (std::size_t i = 0; i < space_dimension; ++i)
711 for (std::size_t j = 0; j < value_size; ++j)
713 u[p * ushape[1] + (j * bs_element + k)]
714 += coefficients[bs_element * i + k] * basis_values(i, j);