365 const std::vector<std::reference_wrapper<const ufcx_form>>& ufcx_forms,
367 const std::vector<std::shared_ptr<
const Function<T, U>>>& coefficients,
368 const std::vector<std::shared_ptr<
const Constant<T>>>& constants,
371 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
374 std::span<const std::int32_t>>& entity_maps,
377 for (
const ufcx_form& ufcx_form : ufcx_forms)
379 if (ufcx_form.rank != (
int)spaces.size())
380 throw std::runtime_error(
"Wrong number of argument spaces for Form.");
381 if (ufcx_form.num_coefficients != (
int)coefficients.size())
383 throw std::runtime_error(
"Mismatch between number of expected and "
384 "provided Form coefficients.");
386 if (ufcx_form.num_constants != (
int)constants.size())
388 throw std::runtime_error(
389 "Mismatch between number of expected and provided Form constants.");
394 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
396 for (std::size_t i = 0; i < spaces.size(); ++i)
398 assert(spaces[i]->elements(form_idx));
399 if (
auto element_hash
400 = ufcx_forms[form_idx].get().finite_element_hashes[i];
403 != spaces[i]->elements(form_idx)->basix_element().hash())
405 throw std::runtime_error(
406 "Cannot create form. Elements are different to "
407 "those used to compile the form.");
413 if (!
mesh and !spaces.empty())
414 mesh = spaces.front()->mesh();
416 throw std::runtime_error(
"No mesh could be associated with the Form.");
417 for (
auto& V : spaces)
419 if (
mesh != V->mesh() and entity_maps.find(V->mesh()) == entity_maps.end())
421 throw std::runtime_error(
422 "Incompatible mesh. entity_maps must be provided.");
426 auto topology =
mesh->topology();
428 const int tdim = topology->dim();
433 const int* integral_offsets = ufcx_forms[0].get().form_integral_offsets;
434 std::vector<int> num_integrals_type(3);
435 for (
int i = 0; i < 3; ++i)
436 num_integrals_type[i] = integral_offsets[i + 1] - integral_offsets[i];
442 mesh->topology_mutable()->create_entities(tdim - 1);
443 mesh->topology_mutable()->create_connectivity(tdim - 1, tdim);
444 mesh->topology_mutable()->create_connectivity(tdim, tdim - 1);
449 using kern_t = std::function<void(T*,
const T*,
const T*,
const U*,
450 const int*,
const std::uint8_t*,
void*)>;
453 auto check_geometry_hash
454 = [&geo =
mesh->geometry()](
const ufcx_integral& integral,
455 std::size_t cell_idx)
457 if (integral.coordinate_element_hash != geo.cmaps().at(cell_idx).hash())
459 throw std::runtime_error(
460 "Generated integral geometry element does not match mesh geometry: "
461 + std::to_string(integral.coordinate_element_hash) +
", "
462 + std::to_string(geo.cmaps().at(cell_idx).hash()));
467 bool needs_facet_permutations =
false;
469 std::vector<std::int32_t> default_cells;
470 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
471 + integral_offsets[
cell],
472 num_integrals_type[
cell]);
474 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
476 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
477 for (
int i = 0; i < num_integrals_type[
cell]; ++i)
479 const int id = ids[i];
480 ufcx_integral* integral
481 = ufcx_form.form_integrals[integral_offsets[
cell] + i];
483 check_geometry_hash(*integral, form_idx);
486 std::vector<int> active_coeffs;
487 for (
int j = 0; j < ufcx_form.num_coefficients; ++j)
489 if (integral->enabled_coefficients[j])
490 active_coeffs.push_back(j);
494 if constexpr (std::is_same_v<T, float>)
495 k = integral->tabulate_tensor_float32;
496#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
497 else if constexpr (std::is_same_v<T, std::complex<float>>)
499 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
500 const scalar_value_t<T>*,
const int*,
501 const unsigned char*,
void*)
>(
502 integral->tabulate_tensor_complex64);
505 else if constexpr (std::is_same_v<T, double>)
506 k = integral->tabulate_tensor_float64;
507#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
508 else if constexpr (std::is_same_v<T, std::complex<double>>)
510 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
511 const scalar_value_t<T>*,
const int*,
512 const unsigned char*,
void*)
>(
513 integral->tabulate_tensor_complex128);
519 throw std::runtime_error(
520 "UFCx kernel function is NULL. Check requested types.");
527 assert(topology->index_maps(tdim).at(form_idx));
528 default_cells.resize(
529 topology->index_maps(tdim).at(form_idx)->size_local(), 0);
530 std::iota(default_cells.begin(), default_cells.end(), 0);
532 {k, default_cells, active_coeffs}});
534 else if (sd != subdomains.end())
537 auto it = std::ranges::lower_bound(sd->second,
id, std::less<>{},
538 [](
auto& a) { return a.first; });
539 if (it != sd->second.end() and it->first ==
id)
543 std::vector<std::int32_t>(it->second.begin(),
549 if (integral->needs_facet_permutations)
550 needs_facet_permutations =
true;
556 std::vector<std::int32_t> default_facets_ext;
558 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
562 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
564 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
567 const int id = ids[i];
568 ufcx_integral* integral
571 check_geometry_hash(*integral, form_idx);
573 std::vector<int> active_coeffs;
574 for (
int j = 0; j < ufcx_form.num_coefficients; ++j)
576 if (integral->enabled_coefficients[j])
577 active_coeffs.push_back(j);
581 if constexpr (std::is_same_v<T, float>)
582 k = integral->tabulate_tensor_float32;
583#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
584 else if constexpr (std::is_same_v<T, std::complex<float>>)
586 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
587 const scalar_value_t<T>*,
const int*,
588 const unsigned char*,
void*)
>(
589 integral->tabulate_tensor_complex64);
592 else if constexpr (std::is_same_v<T, double>)
593 k = integral->tabulate_tensor_float64;
594#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
595 else if constexpr (std::is_same_v<T, std::complex<double>>)
597 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
598 const scalar_value_t<T>*,
const int*,
599 const unsigned char*,
void*)
>(
600 integral->tabulate_tensor_complex128);
607 auto f_to_c = topology->connectivity(tdim - 1, tdim);
609 auto c_to_f = topology->connectivity(tdim, tdim - 1);
614 default_facets_ext.reserve(2 * bfacets.size());
615 for (std::int32_t f : bfacets)
620 default_facets_ext.insert(default_facets_ext.end(), pair.begin(),
624 {k, default_facets_ext, active_coeffs}});
626 else if (sd != subdomains.end())
629 auto it = std::ranges::lower_bound(sd->second,
id, std::less<>{},
630 [](
auto& a) { return a.first; });
631 if (it != sd->second.end() and it->first ==
id)
635 std::vector<std::int32_t>(it->second.begin(),
641 if (integral->needs_facet_permutations)
642 needs_facet_permutations =
true;
648 std::vector<std::int32_t> default_facets_int;
650 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
654 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
656 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
659 std::vector<std::int8_t> interprocess_marker;
662 assert(topology->index_map(tdim - 1));
663 const std::vector<std::int32_t>& interprocess_facets
664 = topology->interprocess_facets();
665 std::int32_t num_facets = topology->index_map(tdim - 1)->size_local()
666 + topology->index_map(tdim - 1)->num_ghosts();
667 interprocess_marker.resize(num_facets, 0);
668 std::ranges::for_each(interprocess_facets,
669 [&interprocess_marker](
auto f)
670 { interprocess_marker[f] = 1; });
675 const int id = ids[i];
676 ufcx_integral* integral
679 check_geometry_hash(*integral, form_idx);
681 std::vector<int> active_coeffs;
682 for (
int j = 0; j < ufcx_form.num_coefficients; ++j)
684 if (integral->enabled_coefficients[j])
685 active_coeffs.push_back(j);
689 if constexpr (std::is_same_v<T, float>)
690 k = integral->tabulate_tensor_float32;
691#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
692 else if constexpr (std::is_same_v<T, std::complex<float>>)
694 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
695 const scalar_value_t<T>*,
const int*,
696 const unsigned char*,
void*)
>(
697 integral->tabulate_tensor_complex64);
700 else if constexpr (std::is_same_v<T, double>)
701 k = integral->tabulate_tensor_float64;
702#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
703 else if constexpr (std::is_same_v<T, std::complex<double>>)
705 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
706 const scalar_value_t<T>*,
const int*,
707 const unsigned char*,
void*)
>(
708 integral->tabulate_tensor_complex128);
714 auto f_to_c = topology->connectivity(tdim - 1, tdim);
716 auto c_to_f = topology->connectivity(tdim, tdim - 1);
721 assert(topology->index_map(tdim - 1));
722 std::int32_t num_facets = topology->index_map(tdim - 1)->size_local();
723 default_facets_int.reserve(4 * num_facets);
724 for (std::int32_t f = 0; f < num_facets; ++f)
726 if (f_to_c->num_links(f) == 2)
730 default_facets_int.insert(default_facets_int.end(), pairs.begin(),
733 else if (interprocess_marker[f])
735 throw std::runtime_error(
736 "Cannot compute interior facet integral over interprocess "
737 "facet. Please use ghost mode shared facet when creating the "
742 {k, default_facets_int, active_coeffs}});
744 else if (sd != subdomains.end())
746 auto it = std::ranges::lower_bound(sd->second,
id, std::less{},
747 [](
auto& a) { return a.first; });
748 if (it != sd->second.end() and it->first ==
id)
752 std::vector<std::int32_t>(it->second.begin(),
758 if (integral->needs_facet_permutations)
759 needs_facet_permutations =
true;
764 return Form<T, U>(spaces, std::move(integrals),
mesh, coefficients, constants,
765 needs_facet_permutations, entity_maps);