188 throw std::runtime_error(
189 "Cannot create sparsity pattern. Form is not a bilinear.");
192 std::shared_ptr
mesh = a.mesh();
194 std::shared_ptr mesh0 = a.function_spaces().at(0)->mesh();
196 std::shared_ptr mesh1 = a.function_spaces().at(1)->mesh();
199 const std::set<IntegralType> types = a.integral_types();
204 int tdim =
mesh->topology()->dim();
205 mesh->topology_mutable()->create_entities(tdim - 1);
206 mesh->topology_mutable()->create_connectivity(tdim - 1, tdim);
211 auto extract_cells = [](std::span<const std::int32_t> facets)
213 assert(facets.size() % 2 == 0);
214 std::vector<std::int32_t> cells;
215 cells.reserve(facets.size() / 2);
216 for (std::size_t i = 0; i < facets.size(); i += 2)
217 cells.push_back(facets[i]);
221 const int num_cell_types =
mesh->topology()->cell_types().size();
222 for (
int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
224 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
225 *a.function_spaces().at(0)->dofmaps(cell_type_idx),
226 *a.function_spaces().at(1)->dofmaps(cell_type_idx)};
229 for (
auto type : types)
234 for (
int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
237 {a.domain_arg(type, 0, i, cell_type_idx),
238 a.domain_arg(type, 1, i, cell_type_idx)},
239 {{dofmaps[0], dofmaps[1]}});
243 for (
int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
247 {extract_cells(a.domain_arg(type, 0, i, 0)),
248 extract_cells(a.domain_arg(type, 1, i, 0))},
249 {{dofmaps[0], dofmaps[1]}});
253 for (
int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
256 {extract_cells(a.domain_arg(type, 0, i, 0)),
257 extract_cells(a.domain_arg(type, 1, i, 0))},
258 {{dofmaps[0], dofmaps[1]}});
262 throw std::runtime_error(
"Unsupported integral type");
369 const std::vector<std::reference_wrapper<const ufcx_form>>& ufcx_forms,
371 const std::vector<std::shared_ptr<
const Function<T, U>>>& coefficients,
372 const std::vector<std::shared_ptr<
const Constant<T>>>& constants,
375 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
377 const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
381 for (
const ufcx_form& ufcx_form : ufcx_forms)
383 if (ufcx_form.rank != (
int)spaces.size())
384 throw std::runtime_error(
"Wrong number of argument spaces for Form.");
385 if (ufcx_form.num_coefficients != (
int)coefficients.size())
387 throw std::runtime_error(
"Mismatch between number of expected and "
388 "provided Form coefficients.");
392 if (ufcx_form.num_constants != (
int)constants.size())
394 throw std::runtime_error(std::format(
395 "Mismatch between number of expected and "
396 "provided Form Constants. Expected {} constants, but got {}.",
397 ufcx_form.num_constants, constants.size()));
399 for (std::size_t c = 0; c < constants.size(); ++c)
401 if (ufcx_form.constant_ranks[c] != (
int)constants[c]->shape.size())
403 throw std::runtime_error(std::format(
404 "Mismatch between expected and actual rank of "
405 "Form Constant. Rank of Constant {} should be {}, but got rank {}.",
406 c, ufcx_form.constant_ranks[c], constants[c]->shape.size()));
408 if (!std::equal(constants[c]->shape.begin(), constants[c]->shape.end(),
409 ufcx_form.constant_shapes[c]))
411 throw std::runtime_error(
412 std::format(
"Mismatch between expected and actual shape of Form "
413 "Constant for Constant {}.",
420 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
422 for (std::size_t i = 0; i < spaces.size(); ++i)
424 assert(spaces[i]->elements(form_idx));
425 if (
auto element_hash
426 = ufcx_forms[form_idx].get().finite_element_hashes[i];
429 != spaces[i]->elements(form_idx)->basix_element().hash())
431 throw std::runtime_error(
432 "Cannot create form. Elements are different to "
433 "those used to compile the form.");
439 if (!
mesh and !spaces.empty())
440 mesh = spaces.front()->mesh();
442 throw std::runtime_error(
"No mesh could be associated with the Form.");
444 auto topology =
mesh->topology();
446 const int tdim = topology->dim();
451 const int* integral_offsets = ufcx_forms[0].get().form_integral_offsets;
452 std::vector<int> num_integrals_type(3);
453 for (
int i = 0; i < 3; ++i)
454 num_integrals_type[i] = integral_offsets[i + 1] - integral_offsets[i];
460 mesh->topology_mutable()->create_entities(tdim - 1);
461 mesh->topology_mutable()->create_connectivity(tdim - 1, tdim);
462 mesh->topology_mutable()->create_connectivity(tdim, tdim - 1);
467 using kern_t = std::function<void(T*,
const T*,
const T*,
const U*,
468 const int*,
const std::uint8_t*,
void*)>;
471 auto check_geometry_hash
472 = [&geo =
mesh->geometry()](
const ufcx_integral& integral,
473 std::size_t cell_idx)
475 if (integral.coordinate_element_hash != geo.cmaps().at(cell_idx).hash())
477 throw std::runtime_error(
478 "Generated integral geometry element does not match mesh geometry: "
479 + std::to_string(integral.coordinate_element_hash) +
", "
480 + std::to_string(geo.cmaps().at(cell_idx).hash()));
485 bool needs_facet_permutations =
false;
487 std::vector<std::int32_t> default_cells;
488 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
489 + integral_offsets[
cell],
490 num_integrals_type[
cell]);
492 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
494 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
495 for (
int i = 0; i < num_integrals_type[
cell]; ++i)
497 const int id = ids[i];
498 ufcx_integral* integral
499 = ufcx_form.form_integrals[integral_offsets[
cell] + i];
501 check_geometry_hash(*integral, form_idx);
504 std::vector<int> active_coeffs;
505 for (
int j = 0; j < ufcx_form.num_coefficients; ++j)
507 if (integral->enabled_coefficients[j])
508 active_coeffs.push_back(j);
512 if constexpr (std::is_same_v<T, float>)
513 k = integral->tabulate_tensor_float32;
514#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
515 else if constexpr (std::is_same_v<T, std::complex<float>>)
517 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
518 const scalar_value_t<T>*,
const int*,
519 const unsigned char*,
void*)
>(
520 integral->tabulate_tensor_complex64);
523 else if constexpr (std::is_same_v<T, double>)
524 k = integral->tabulate_tensor_float64;
525#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
526 else if constexpr (std::is_same_v<T, std::complex<double>>)
528 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
529 const scalar_value_t<T>*,
const int*,
530 const unsigned char*,
void*)
>(
531 integral->tabulate_tensor_complex128);
537 throw std::runtime_error(
538 "UFCx kernel function is NULL. Check requested types.");
545 assert(topology->index_maps(tdim).at(form_idx));
546 default_cells.resize(
547 topology->index_maps(tdim).at(form_idx)->size_local(), 0);
548 std::iota(default_cells.begin(), default_cells.end(), 0);
550 {k, default_cells, active_coeffs}});
552 else if (sd != subdomains.end())
555 auto it = std::ranges::lower_bound(sd->second,
id, std::less<>{},
556 [](
auto& a) { return a.first; });
557 if (it != sd->second.end() and it->first ==
id)
561 std::vector<std::int32_t>(it->second.begin(),
567 if (integral->needs_facet_permutations)
568 needs_facet_permutations =
true;
574 std::vector<std::int32_t> default_facets_ext;
576 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
580 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
582 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
585 const int id = ids[i];
586 ufcx_integral* integral
589 check_geometry_hash(*integral, form_idx);
591 std::vector<int> active_coeffs;
592 for (
int j = 0; j < ufcx_form.num_coefficients; ++j)
594 if (integral->enabled_coefficients[j])
595 active_coeffs.push_back(j);
599 if constexpr (std::is_same_v<T, float>)
600 k = integral->tabulate_tensor_float32;
601#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
602 else if constexpr (std::is_same_v<T, std::complex<float>>)
604 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
605 const scalar_value_t<T>*,
const int*,
606 const unsigned char*,
void*)
>(
607 integral->tabulate_tensor_complex64);
610 else if constexpr (std::is_same_v<T, double>)
611 k = integral->tabulate_tensor_float64;
612#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
613 else if constexpr (std::is_same_v<T, std::complex<double>>)
615 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
616 const scalar_value_t<T>*,
const int*,
617 const unsigned char*,
void*)
>(
618 integral->tabulate_tensor_complex128);
625 auto f_to_c = topology->connectivity(tdim - 1, tdim);
627 auto c_to_f = topology->connectivity(tdim, tdim - 1);
632 default_facets_ext.reserve(2 * bfacets.size());
633 for (std::int32_t f : bfacets)
636 std::array<std::int32_t, 2> pair
638 default_facets_ext.insert(default_facets_ext.end(), pair.begin(),
642 {k, default_facets_ext, active_coeffs}});
644 else if (sd != subdomains.end())
647 auto it = std::ranges::lower_bound(sd->second,
id, std::less<>{},
648 [](
auto& a) { return a.first; });
649 if (it != sd->second.end() and it->first ==
id)
653 std::vector<std::int32_t>(it->second.begin(),
659 if (integral->needs_facet_permutations)
660 needs_facet_permutations =
true;
666 std::vector<std::int32_t> default_facets_int;
668 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
672 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
674 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
677 std::vector<std::int8_t> interprocess_marker;
680 assert(topology->index_map(tdim - 1));
681 const std::vector<std::int32_t>& interprocess_facets
682 = topology->interprocess_facets();
683 std::int32_t num_facets = topology->index_map(tdim - 1)->size_local()
684 + topology->index_map(tdim - 1)->num_ghosts();
685 interprocess_marker.resize(num_facets, 0);
686 std::ranges::for_each(interprocess_facets,
687 [&interprocess_marker](
auto f)
688 { interprocess_marker[f] = 1; });
693 const int id = ids[i];
694 ufcx_integral* integral
697 check_geometry_hash(*integral, form_idx);
699 std::vector<int> active_coeffs;
700 for (
int j = 0; j < ufcx_form.num_coefficients; ++j)
702 if (integral->enabled_coefficients[j])
703 active_coeffs.push_back(j);
707 if constexpr (std::is_same_v<T, float>)
708 k = integral->tabulate_tensor_float32;
709#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
710 else if constexpr (std::is_same_v<T, std::complex<float>>)
712 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
713 const scalar_value_t<T>*,
const int*,
714 const unsigned char*,
void*)
>(
715 integral->tabulate_tensor_complex64);
718 else if constexpr (std::is_same_v<T, double>)
719 k = integral->tabulate_tensor_float64;
720#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
721 else if constexpr (std::is_same_v<T, std::complex<double>>)
723 k =
reinterpret_cast<void (*)(T*,
const T*,
const T*,
724 const scalar_value_t<T>*,
const int*,
725 const unsigned char*,
void*)
>(
726 integral->tabulate_tensor_complex128);
732 auto f_to_c = topology->connectivity(tdim - 1, tdim);
734 auto c_to_f = topology->connectivity(tdim, tdim - 1);
739 assert(topology->index_map(tdim - 1));
740 std::int32_t num_facets = topology->index_map(tdim - 1)->size_local();
741 default_facets_int.reserve(4 * num_facets);
742 for (std::int32_t f = 0; f < num_facets; ++f)
744 if (f_to_c->num_links(f) == 2)
746 std::array<std::int32_t, 4> pairs
748 default_facets_int.insert(default_facets_int.end(), pairs.begin(),
751 else if (interprocess_marker[f])
753 throw std::runtime_error(
754 "Cannot compute interior facet integral over interprocess "
755 "facet. Please use ghost mode shared facet when creating the "
760 {k, default_facets_int, active_coeffs}});
762 else if (sd != subdomains.end())
764 auto it = std::ranges::lower_bound(sd->second,
id, std::less{},
765 [](
auto& a) { return a.first; });
766 if (it != sd->second.end() and it->first ==
id)
770 std::vector<std::int32_t>(it->second.begin(),
776 if (integral->needs_facet_permutations)
777 needs_facet_permutations =
true;
782 return Form<T, U>(spaces, std::move(integrals),
mesh, coefficients, constants,
783 needs_facet_permutations, entity_maps);