DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
utils.h
Go to the documentation of this file.
1// Copyright (C) 2013-2020 Johan Hake, Jan Blechta and Garth N. Wells
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "Constant.h"
10#include "CoordinateElement.h"
11#include "DofMap.h"
12#include "ElementDofLayout.h"
13#include "Expression.h"
14#include "Form.h"
15#include "Function.h"
16#include "FunctionSpace.h"
17#include "sparsitybuild.h"
18#include <algorithm>
19#include <array>
20#include <concepts>
21#include <dolfinx/common/types.h>
22#include <dolfinx/la/SparsityPattern.h>
23#include <dolfinx/mesh/EntityMap.h>
24#include <dolfinx/mesh/Topology.h>
25#include <dolfinx/mesh/cell_types.h>
26#include <dolfinx/mesh/utils.h>
27#include <format>
28#include <functional>
29#include <memory>
30#include <optional>
31#include <set>
32#include <span>
33#include <stdexcept>
34#include <string>
35#include <tuple>
36#include <type_traits>
37#include <ufcx.h>
38#include <utility>
39#include <vector>
40
43namespace basix
44{
45template <std::floating_point T>
46class FiniteElement;
47}
48
49namespace dolfinx::common
50{
51class IndexMap;
52}
53
54namespace dolfinx::fem
55{
56template <dolfinx::scalar T, std::floating_point U>
57class Expression;
58
59namespace impl
60{
67template <int num_cells>
68std::array<std::int32_t, 2 * num_cells>
69get_cell_facet_pairs(std::int32_t f, std::span<const std::int32_t> cells,
71{
72 // Loop over cells sharing facet
73 assert(cells.size() == num_cells);
74 std::array<std::int32_t, 2 * num_cells> cell_local_facet_pairs;
75 for (int c = 0; c < num_cells; ++c)
76 {
77 // Get local index of facet with respect to the cell
78 std::int32_t cell = cells[c];
79 auto cell_facets = c_to_f.links(cell);
80 auto facet_it = std::find(cell_facets.begin(), cell_facets.end(), f);
81 assert(facet_it != cell_facets.end());
82 int local_f = std::distance(cell_facets.begin(), facet_it);
83 cell_local_facet_pairs[2 * c] = cell;
84 cell_local_facet_pairs[2 * c + 1] = local_f;
85 }
86
87 return cell_local_facet_pairs;
88}
89} // namespace impl
90
117std::vector<std::int32_t>
119 const mesh::Topology& topology,
120 std::span<const std::int32_t> entities);
121
129template <dolfinx::scalar T, std::floating_point U>
130std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>>
131extract_function_spaces(const std::vector<std::vector<const Form<T, U>*>>& a)
132{
133 std::vector<
134 std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>>
135 spaces(
136 a.size(),
137 std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>(
138 a.front().size()));
139 for (std::size_t i = 0; i < a.size(); ++i)
140 {
141 for (std::size_t j = 0; j < a[i].size(); ++j)
142 {
143 if (const Form<T, U>* form = a[i][j]; form)
144 spaces[i][j] = {form->function_spaces()[0], form->function_spaces()[1]};
145 }
146 }
147 return spaces;
148}
149
155template <dolfinx::scalar T, std::floating_point U>
157{
158 std::shared_ptr mesh = a.mesh();
159 assert(mesh);
160
161 // Get index maps and block sizes from the DOF maps. Note that in
162 // mixed-topology meshes, despite there being multiple DOF maps, the
163 // index maps and block sizes are the same.
164 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
165 *a.function_spaces().at(0)->dofmaps(0),
166 *a.function_spaces().at(1)->dofmaps(0)};
167
168 const std::array index_maps{dofmaps[0].get().index_map,
169 dofmaps[1].get().index_map};
170 const std::array bs
171 = {dofmaps[0].get().index_map_bs(), dofmaps[1].get().index_map_bs()};
172
173 la::SparsityPattern pattern(mesh->comm(), index_maps, bs);
174 build_sparsity_pattern(pattern, a);
175 return pattern;
176}
177
183template <dolfinx::scalar T, std::floating_point U>
185{
186 if (a.rank() != 2)
187 {
188 throw std::runtime_error(
189 "Cannot create sparsity pattern. Form is not a bilinear.");
190 }
191
192 std::shared_ptr mesh = a.mesh();
193 assert(mesh);
194 std::shared_ptr mesh0 = a.function_spaces().at(0)->mesh();
195 assert(mesh0);
196 std::shared_ptr mesh1 = a.function_spaces().at(1)->mesh();
197 assert(mesh1);
198
199 const std::set<IntegralType> types = a.integral_types();
200 if (types.find(IntegralType::interior_facet) != types.end()
201 or types.find(IntegralType::exterior_facet) != types.end())
202 {
203 // FIXME: cleanup these calls? Some of the happen internally again.
204 int tdim = mesh->topology()->dim();
205 mesh->topology_mutable()->create_entities(tdim - 1);
206 mesh->topology_mutable()->create_connectivity(tdim - 1, tdim);
207 }
208
209 common::Timer t0("Build sparsity");
210
211 auto extract_cells = [](std::span<const std::int32_t> facets)
212 {
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]);
218 return cells;
219 };
220
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)
223 {
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)};
227
228 // Create and build sparsity pattern
229 for (auto type : types)
230 {
231 switch (type)
232 {
234 for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
235 {
236 sparsitybuild::cells(pattern,
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]}});
240 }
241 break;
243 for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
244 {
246 pattern,
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]}});
250 }
251 break;
253 for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i)
254 {
255 sparsitybuild::cells(pattern,
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]}});
259 }
260 break;
261 default:
262 throw std::runtime_error("Unsupported integral type");
263 }
264 }
265 }
266
267 t0.stop();
268}
269
271template <std::floating_point T>
273 const std::vector<int>& parent_map
274 = {})
275{
276 // Create subdofmaps and compute offset
277 std::vector<int> offsets(1, 0);
278 std::vector<dolfinx::fem::ElementDofLayout> sub_doflayout;
279 int bs = element.block_size();
280 for (int i = 0; i < element.num_sub_elements(); ++i)
281 {
282 // The ith sub-element. For mixed elements this is subelements()[i]. For
283 // blocked elements, the sub-element will always be the same, so we'll use
284 // sub_elements()[0]
285 std::shared_ptr<const fem::FiniteElement<T>> sub_e
286 = element.sub_elements()[bs > 1 ? 0 : i];
287
288 // In a mixed element DOFs are ordered element by element, so the offset to
289 // the next sub-element is sub_e->space_dimension(). Blocked elements use
290 // xxyyzz ordering, so the offset to the next sub-element is 1
291
292 std::vector<int> parent_map_sub(sub_e->space_dimension(), offsets.back());
293 for (std::size_t j = 0; j < parent_map_sub.size(); ++j)
294 parent_map_sub[j] += bs * j;
295 offsets.push_back(offsets.back() + (bs > 1 ? 1 : sub_e->space_dimension()));
296 sub_doflayout.push_back(
297 dolfinx::fem::create_element_dof_layout(*sub_e, parent_map_sub));
298 }
299
300 return ElementDofLayout(bs, element.entity_dofs(),
301 element.entity_closure_dofs(), parent_map,
302 sub_doflayout);
303}
304
313DofMap
314create_dofmap(MPI_Comm comm, const ElementDofLayout& layout,
315 mesh::Topology& topology,
316 const std::function<void(std::span<std::int32_t>, std::uint32_t)>&
317 permute_inv,
318 const std::function<std::vector<int>(
319 const graph::AdjacencyList<std::int32_t>&)>& reorder_fn);
320
331std::vector<DofMap> create_dofmaps(
332 MPI_Comm comm, const std::vector<ElementDofLayout>& layouts,
333 mesh::Topology& topology,
334 const std::function<void(std::span<std::int32_t>, std::uint32_t)>&
335 permute_inv,
336 const std::function<std::vector<int>(
337 const graph::AdjacencyList<std::int32_t>&)>& reorder_fn);
338
342std::vector<std::string> get_coefficient_names(const ufcx_form& ufcx_form);
343
347std::vector<std::string> get_constant_names(const ufcx_form& ufcx_form);
348
367template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
369 const std::vector<std::reference_wrapper<const ufcx_form>>& ufcx_forms,
370 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
371 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients,
372 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
373 const std::map<
375 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
376 subdomains,
377 const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
378 entity_maps,
379 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
380{
381 for (const ufcx_form& ufcx_form : ufcx_forms)
382 {
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())
386 {
387 throw std::runtime_error("Mismatch between number of expected and "
388 "provided Form coefficients.");
389 }
390
391 // Check Constants for rank and size consistency
392 if (ufcx_form.num_constants != (int)constants.size())
393 {
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()));
398 }
399 for (std::size_t c = 0; c < constants.size(); ++c)
400 {
401 if (ufcx_form.constant_ranks[c] != (int)constants[c]->shape.size())
402 {
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()));
407 }
408 if (!std::equal(constants[c]->shape.begin(), constants[c]->shape.end(),
409 ufcx_form.constant_shapes[c]))
410 {
411 throw std::runtime_error(
412 std::format("Mismatch between expected and actual shape of Form "
413 "Constant for Constant {}.",
414 c));
415 }
416 }
417 }
418
419 // Check argument function spaces
420 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
421 {
422 for (std::size_t i = 0; i < spaces.size(); ++i)
423 {
424 assert(spaces[i]->elements(form_idx));
425 if (auto element_hash
426 = ufcx_forms[form_idx].get().finite_element_hashes[i];
427 element_hash != 0
428 and element_hash
429 != spaces[i]->elements(form_idx)->basix_element().hash())
430 {
431 throw std::runtime_error(
432 "Cannot create form. Elements are different to "
433 "those used to compile the form.");
434 }
435 }
436 }
437
438 // Extract mesh from FunctionSpace, and check they are the same
439 if (!mesh and !spaces.empty())
440 mesh = spaces.front()->mesh();
441 if (!mesh)
442 throw std::runtime_error("No mesh could be associated with the Form.");
443
444 auto topology = mesh->topology();
445 assert(topology);
446 const int tdim = topology->dim();
447
448 // NOTE: This assumes all forms in mixed-topology meshes have the same
449 // integral offsets. Since the UFL forms for each type of cell should be
450 // the same, I think this assumption is OK.
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];
455
456 // Create facets, if required
457 if (num_integrals_type[exterior_facet] > 0
458 or num_integrals_type[interior_facet] > 0)
459 {
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);
463 }
464
465 // Get list of integral IDs, and load tabulate tensor into memory for
466 // each
467 using kern_t = std::function<void(T*, const T*, const T*, const U*,
468 const int*, const std::uint8_t*, void*)>;
469 std::map<std::tuple<IntegralType, int, int>, integral_data<T, U>> integrals;
470
471 auto check_geometry_hash
472 = [&geo = mesh->geometry()](const ufcx_integral& integral,
473 std::size_t cell_idx)
474 {
475 if (integral.coordinate_element_hash != geo.cmaps().at(cell_idx).hash())
476 {
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()));
481 }
482 };
483
484 // Attach cell kernels
485 bool needs_facet_permutations = false;
486 {
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]);
491 auto sd = subdomains.find(IntegralType::cell);
492 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
493 {
494 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
495 for (int i = 0; i < num_integrals_type[cell]; ++i)
496 {
497 const int id = ids[i];
498 ufcx_integral* integral
499 = ufcx_form.form_integrals[integral_offsets[cell] + i];
500 assert(integral);
501 check_geometry_hash(*integral, form_idx);
502
503 // Build list of active coefficients
504 std::vector<int> active_coeffs;
505 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
506 {
507 if (integral->enabled_coefficients[j])
508 active_coeffs.push_back(j);
509 }
510
511 kern_t k = nullptr;
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>>)
516 {
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);
521 }
522#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
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>>)
527 {
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);
532 }
533#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
534
535 if (!k)
536 {
537 throw std::runtime_error(
538 "UFCx kernel function is NULL. Check requested types.");
539 }
540
541 // Build list of entities to assemble over
542 if (id == -1)
543 {
544 // Default kernel, operates on all (owned) cells
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);
549 integrals.insert({{IntegralType::cell, i, form_idx},
550 {k, default_cells, active_coeffs}});
551 }
552 else if (sd != subdomains.end())
553 {
554 // NOTE: This requires that pairs are sorted
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)
558 {
559 integrals.insert({{IntegralType::cell, i, form_idx},
560 {k,
561 std::vector<std::int32_t>(it->second.begin(),
562 it->second.end()),
563 active_coeffs}});
564 }
565 }
566
567 if (integral->needs_facet_permutations)
568 needs_facet_permutations = true;
569 }
570 }
571 }
572
573 // Attach exterior facet kernels
574 std::vector<std::int32_t> default_facets_ext;
575 {
576 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
577 + integral_offsets[exterior_facet],
578 num_integrals_type[exterior_facet]);
579 auto sd = subdomains.find(IntegralType::exterior_facet);
580 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
581 {
582 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
583 for (int i = 0; i < num_integrals_type[exterior_facet]; ++i)
584 {
585 const int id = ids[i];
586 ufcx_integral* integral
587 = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i];
588 assert(integral);
589 check_geometry_hash(*integral, form_idx);
590
591 std::vector<int> active_coeffs;
592 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
593 {
594 if (integral->enabled_coefficients[j])
595 active_coeffs.push_back(j);
596 }
597
598 kern_t k = nullptr;
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>>)
603 {
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);
608 }
609#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
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>>)
614 {
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);
619 }
620#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
621 assert(k);
622
623 // Build list of entities to assembler over
624 const std::vector bfacets = mesh::exterior_facet_indices(*topology);
625 auto f_to_c = topology->connectivity(tdim - 1, tdim);
626 assert(f_to_c);
627 auto c_to_f = topology->connectivity(tdim, tdim - 1);
628 assert(c_to_f);
629 if (id == -1)
630 {
631 // Default kernel, operates on all (owned) exterior facets
632 default_facets_ext.reserve(2 * bfacets.size());
633 for (std::int32_t f : bfacets)
634 {
635 // There will only be one pair for an exterior facet integral
636 std::array<std::int32_t, 2> pair
637 = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f);
638 default_facets_ext.insert(default_facets_ext.end(), pair.begin(),
639 pair.end());
640 }
641 integrals.insert({{IntegralType::exterior_facet, i, form_idx},
642 {k, default_facets_ext, active_coeffs}});
643 }
644 else if (sd != subdomains.end())
645 {
646 // NOTE: This requires that pairs are sorted
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)
650 {
651 integrals.insert({{IntegralType::exterior_facet, i, form_idx},
652 {k,
653 std::vector<std::int32_t>(it->second.begin(),
654 it->second.end()),
655 active_coeffs}});
656 }
657 }
658
659 if (integral->needs_facet_permutations)
660 needs_facet_permutations = true;
661 }
662 }
663 }
664
665 // Attach interior facet kernels
666 std::vector<std::int32_t> default_facets_int;
667 {
668 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
669 + integral_offsets[interior_facet],
670 num_integrals_type[interior_facet]);
671 auto sd = subdomains.find(IntegralType::interior_facet);
672 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
673 {
674 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
675
676 // Create indicator for interprocess facets
677 std::vector<std::int8_t> interprocess_marker;
678 if (num_integrals_type[interior_facet] > 0)
679 {
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; });
689 }
690
691 for (int i = 0; i < num_integrals_type[interior_facet]; ++i)
692 {
693 const int id = ids[i];
694 ufcx_integral* integral
695 = ufcx_form.form_integrals[integral_offsets[interior_facet] + i];
696 assert(integral);
697 check_geometry_hash(*integral, form_idx);
698
699 std::vector<int> active_coeffs;
700 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
701 {
702 if (integral->enabled_coefficients[j])
703 active_coeffs.push_back(j);
704 }
705
706 kern_t k = nullptr;
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>>)
711 {
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);
716 }
717#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
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>>)
722 {
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);
727 }
728#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
729 assert(k);
730
731 // Build list of entities to assembler over
732 auto f_to_c = topology->connectivity(tdim - 1, tdim);
733 assert(f_to_c);
734 auto c_to_f = topology->connectivity(tdim, tdim - 1);
735 assert(c_to_f);
736 if (id == -1)
737 {
738 // Default kernel, operates on all (owned) interior facets
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)
743 {
744 if (f_to_c->num_links(f) == 2)
745 {
746 std::array<std::int32_t, 4> pairs
747 = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f);
748 default_facets_int.insert(default_facets_int.end(), pairs.begin(),
749 pairs.end());
750 }
751 else if (interprocess_marker[f])
752 {
753 throw std::runtime_error(
754 "Cannot compute interior facet integral over interprocess "
755 "facet. Please use ghost mode shared facet when creating the "
756 "mesh");
757 }
758 }
759 integrals.insert({{IntegralType::interior_facet, i, form_idx},
760 {k, default_facets_int, active_coeffs}});
761 }
762 else if (sd != subdomains.end())
763 {
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)
767 {
768 integrals.insert({{IntegralType::interior_facet, i, form_idx},
769 {k,
770 std::vector<std::int32_t>(it->second.begin(),
771 it->second.end()),
772 active_coeffs}});
773 }
774 }
775
776 if (integral->needs_facet_permutations)
777 needs_facet_permutations = true;
778 }
779 }
780 }
781
782 return Form<T, U>(spaces, std::move(integrals), mesh, coefficients, constants,
783 needs_facet_permutations, entity_maps);
784}
785
800template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
802 const ufcx_form& ufcx_form,
803 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
804 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
805 coefficients,
806 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
807 const std::map<
809 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
810 subdomains,
811 const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
812 entity_maps,
813 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
814{
815 // Place coefficients in appropriate order
816 std::vector<std::shared_ptr<const Function<T, U>>> coeff_map;
817 for (const std::string& name : get_coefficient_names(ufcx_form))
818 {
819 if (auto it = coefficients.find(name); it != coefficients.end())
820 coeff_map.push_back(it->second);
821 else
822 {
823 throw std::runtime_error("Form coefficient \"" + name
824 + "\" not provided.");
825 }
826 }
827
828 // Place constants in appropriate order
829 std::vector<std::shared_ptr<const Constant<T>>> const_map;
830 for (const std::string& name : get_constant_names(ufcx_form))
831 {
832 if (auto it = constants.find(name); it != constants.end())
833 const_map.push_back(it->second);
834 else
835 throw std::runtime_error("Form constant \"" + name + "\" not provided.");
836 }
837
838 return create_form_factory({ufcx_form}, spaces, coeff_map, const_map,
839 subdomains, entity_maps, mesh);
840}
841
860template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
862 ufcx_form* (*fptr)(),
863 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
864 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
865 coefficients,
866 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
867 const std::map<
869 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
870 subdomains,
871 const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
872 entity_maps,
873 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
874{
875 ufcx_form* form = fptr();
876 Form<T, U> L = create_form<T, U>(*form, spaces, coefficients, constants,
877 subdomains, entity_maps, mesh);
878 std::free(form);
879 return L;
880}
881
883template <std::floating_point T>
885 std::shared_ptr<mesh::Mesh<T>> mesh,
886 std::shared_ptr<const fem::FiniteElement<T>> e,
887 std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
888 reorder_fn
889 = nullptr)
890{
891 // TODO: check cell type of e (need to add method to fem::FiniteElement)
892 assert(e);
893 assert(mesh);
894 assert(mesh->topology());
895 if (e->cell_type() != mesh->topology()->cell_type())
896 throw std::runtime_error("Cell type of element and mesh must match.");
897
898 // Create element dof layout
900
901 // Create a dofmap
902 std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv
903 = e->needs_dof_permutations() ? e->dof_permutation_fn(true, true)
904 : nullptr;
905 auto dofmap = std::make_shared<const DofMap>(create_dofmap(
906 mesh->comm(), layout, *mesh->topology(), permute_inv, reorder_fn));
907
908 return FunctionSpace(mesh, e, dofmap);
909}
910
912template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
914 const ufcx_expression& e,
915 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients,
916 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
917 std::shared_ptr<const FunctionSpace<U>> argument_space = nullptr)
918{
919 if (!coefficients.empty())
920 {
921 assert(coefficients.front());
922 assert(coefficients.front()->function_space());
923 std::shared_ptr<const mesh::Mesh<U>> mesh
924 = coefficients.front()->function_space()->mesh();
925 if (mesh->geometry().cmap().hash() != e.coordinate_element_hash)
926 {
927 throw std::runtime_error(
928 "Expression and mesh geometric maps do not match.");
929 }
930 }
931
932 if (e.rank > 0 and !argument_space)
933 {
934 throw std::runtime_error("Expression has Argument but no Argument "
935 "function space was provided.");
936 }
937
938 std::vector<U> X(e.points, e.points + e.num_points * e.entity_dimension);
939 std::array<std::size_t, 2> Xshape
940 = {static_cast<std::size_t>(e.num_points),
941 static_cast<std::size_t>(e.entity_dimension)};
942 std::vector<std::size_t> value_shape(e.value_shape,
943 e.value_shape + e.num_components);
944 std::function<void(T*, const T*, const T*, const scalar_value_t<T>*,
945 const int*, const std::uint8_t*, void*)>
946 tabulate_tensor = nullptr;
947 if constexpr (std::is_same_v<T, float>)
948 tabulate_tensor = e.tabulate_tensor_float32;
949#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
950 else if constexpr (std::is_same_v<T, std::complex<float>>)
951 {
952 tabulate_tensor = reinterpret_cast<void (*)(
953 T*, const T*, const T*, const scalar_value_t<T>*, const int*,
954 const unsigned char*, void*)>(e.tabulate_tensor_complex64);
955 }
956#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
957 else if constexpr (std::is_same_v<T, double>)
958 tabulate_tensor = e.tabulate_tensor_float64;
959#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
960 else if constexpr (std::is_same_v<T, std::complex<double>>)
961 {
962 tabulate_tensor = reinterpret_cast<void (*)(
963 T*, const T*, const T*, const scalar_value_t<T>*, const int*,
964 const unsigned char*, void*)>(e.tabulate_tensor_complex128);
965 }
966#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
967 else
968 throw std::runtime_error("Type not supported.");
969
970 assert(tabulate_tensor);
971 return Expression(coefficients, constants, std::span<const U>(X), Xshape,
972 tabulate_tensor, value_shape, argument_space);
973}
974
977template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
979 const ufcx_expression& e,
980 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
981 coefficients,
982 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
983 std::shared_ptr<const FunctionSpace<U>> argument_space = nullptr)
984{
985 // Place coefficients in appropriate order
986 std::vector<std::shared_ptr<const Function<T, U>>> coeff_map;
987 std::vector<std::string> coefficient_names;
988 coefficient_names.reserve(e.num_coefficients);
989 for (int i = 0; i < e.num_coefficients; ++i)
990 coefficient_names.push_back(e.coefficient_names[i]);
991
992 for (const std::string& name : coefficient_names)
993 {
994 if (auto it = coefficients.find(name); it != coefficients.end())
995 coeff_map.push_back(it->second);
996 else
997 {
998 throw std::runtime_error("Expression coefficient \"" + name
999 + "\" not provided.");
1000 }
1001 }
1002
1003 // Place constants in appropriate order
1004 std::vector<std::shared_ptr<const Constant<T>>> const_map;
1005 std::vector<std::string> constant_names;
1006 constant_names.reserve(e.num_constants);
1007 for (int i = 0; i < e.num_constants; ++i)
1008 constant_names.push_back(e.constant_names[i]);
1009
1010 for (const std::string& name : constant_names)
1011 {
1012 if (auto it = constants.find(name); it != constants.end())
1013 const_map.push_back(it->second);
1014 else
1015 {
1016 throw std::runtime_error("Expression constant \"" + name
1017 + "\" not provided.");
1018 }
1019 }
1020
1021 return create_expression(e, coeff_map, const_map, argument_space);
1022}
1023
1024} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
Definition CoordinateElement.h:26
Definition IndexMap.h:97
Timer for measuring and logging elapsed time durations.
Definition Timer.h:40
std::chrono::duration< double, Period > stop()
Stop timer and return elapsed time.
Definition Timer.h:92
Constant (in space) value which can be attached to a Form.
Definition Constant.h:22
Degree-of-freedom map.
Definition DofMap.h:73
Definition ElementDofLayout.h:31
An Expression represents a mathematical expression evaluated at a pre-defined points on a reference c...
Definition Expression.h:41
Model of a finite element.
Definition FiniteElement.h:57
const std::vector< std::shared_ptr< const FiniteElement< geometry_type > > > & sub_elements() const noexcept
Get subelements (if any)
Definition FiniteElement.cpp:396
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const noexcept
Local DOFs associated with each sub-entity of the cell.
Definition FiniteElement.cpp:340
int num_sub_elements() const noexcept
Number of sub elements (for a mixed or blocked element).
Definition FiniteElement.cpp:383
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const noexcept
Local DOFs associated with the closure of each sub-entity of the cell.
Definition FiniteElement.cpp:347
int block_size() const noexcept
Block size of the finite element function space.
Definition FiniteElement.cpp:359
A representation of finite element variational forms.
Definition Form.h:116
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
Definition Function.h:47
This class provides a static adjacency list data structure.
Definition AdjacencyList.h:38
std::span< LinkData > links(std::size_t node)
Get the links (edges) for given node.
Definition AdjacencyList.h:156
Definition SparsityPattern.h:26
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:41
std::array< std::int32_t, 2 *num_cells > get_cell_facet_pairs(std::int32_t f, std::span< const std::int32_t > cells, const graph::AdjacencyList< std::int32_t > &c_to_f)
Definition utils.h:69
Functions supporting mesh operations.
Miscellaneous classes, functions and types.
Definition dolfinx_common.h:8
void interior_facets(la::SparsityPattern &pattern, std::array< std::span< const std::int32_t >, 2 > cells, std::array< std::reference_wrapper< const DofMap >, 2 > dofmaps)
Iterate over interior facets and insert entries into sparsity pattern.
Definition sparsitybuild.cpp:28
void cells(la::SparsityPattern &pattern, std::array< std::span< const std::int32_t >, 2 > cells, std::array< std::reference_wrapper< const DofMap >, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition sparsitybuild.cpp:16
Finite element method functionality.
Definition assemble_expression_impl.h:23
std::vector< std::string > get_constant_names(const ufcx_form &ufcx_form)
Get the name of each constant in a UFC form.
Definition utils.cpp:127
FunctionSpace< T > create_functionspace(std::shared_ptr< mesh::Mesh< T > > mesh, std::shared_ptr< const fem::FiniteElement< T > > e, std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> reorder_fn=nullptr)
NEW Create a function space from a fem::FiniteElement.
Definition utils.h:884
Form< T, U > create_form_factory(const std::vector< std::reference_wrapper< const ufcx_form > > &ufcx_forms, const std::vector< std::shared_ptr< const FunctionSpace< U > > > &spaces, const std::vector< std::shared_ptr< const Function< T, U > > > &coefficients, const std::vector< std::shared_ptr< const Constant< T > > > &constants, const std::map< IntegralType, std::vector< std::pair< std::int32_t, std::span< const std::int32_t > > > > &subdomains, const std::vector< std::reference_wrapper< const mesh::EntityMap > > &entity_maps, std::shared_ptr< const mesh::Mesh< U > > mesh=nullptr)
Create a Form from UFCx input with coefficients and constants passed in the required order.
Definition utils.h:368
std::vector< std::string > get_coefficient_names(const ufcx_form &ufcx_form)
Definition utils.cpp:120
std::vector< std::vector< std::array< std::shared_ptr< const FunctionSpace< U > >, 2 > > > extract_function_spaces(const std::vector< std::vector< const Form< T, U > * > > &a)
Extract test (0) and trial (1) function spaces pairs for each bilinear form for a rectangular array o...
Definition utils.h:131
Expression< T, U > create_expression(const ufcx_expression &e, const std::vector< std::shared_ptr< const Function< T, U > > > &coefficients, const std::vector< std::shared_ptr< const Constant< T > > > &constants, std::shared_ptr< const FunctionSpace< U > > argument_space=nullptr)
Create Expression from UFC.
Definition utils.h:913
std::vector< std::int32_t > compute_integration_domains(IntegralType integral_type, const mesh::Topology &topology, std::span< const std::int32_t > entities)
Given an integral type and a set of entities, computes and return data for the entities that should b...
Definition utils.cpp:135
Form< T, U > create_form(const ufcx_form &ufcx_form, const std::vector< std::shared_ptr< const FunctionSpace< U > > > &spaces, const std::map< std::string, std::shared_ptr< const Function< T, U > > > &coefficients, const std::map< std::string, std::shared_ptr< const Constant< T > > > &constants, const std::map< IntegralType, std::vector< std::pair< std::int32_t, std::span< const std::int32_t > > > > &subdomains, const std::vector< std::reference_wrapper< const mesh::EntityMap > > &entity_maps, std::shared_ptr< const mesh::Mesh< U > > mesh=nullptr)
Create a Form from UFC input with coefficients and constants resolved by name.
Definition utils.h:801
ElementDofLayout create_element_dof_layout(const fem::FiniteElement< T > &element, const std::vector< int > &parent_map={})
Create an ElementDofLayout from a FiniteElement.
Definition utils.h:272
IntegralType
Type of integral.
Definition Form.h:39
@ interior_facet
Interior facet.
Definition Form.h:42
@ cell
Cell.
Definition Form.h:40
@ exterior_facet
Exterior facet.
Definition Form.h:41
la::SparsityPattern create_sparsity_pattern(const Form< T, U > &a)
Create a sparsity pattern for a given form.
Definition utils.h:156
DofMap create_dofmap(MPI_Comm comm, const ElementDofLayout &layout, mesh::Topology &topology, const std::function< void(std::span< std::int32_t >, std::uint32_t)> &permute_inv, const std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> &reorder_fn)
Create a dof map on mesh.
Definition utils.cpp:31
std::vector< DofMap > create_dofmaps(MPI_Comm comm, const std::vector< ElementDofLayout > &layouts, mesh::Topology &topology, const std::function< void(std::span< std::int32_t >, std::uint32_t)> &permute_inv, const std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> &reorder_fn)
Create a set of dofmaps on a given topology.
Definition utils.cpp:69
FunctionSpace(U mesh, V element, W dofmap) -> FunctionSpace< typename std::remove_cvref< typename U::element_type >::type::geometry_type::value_type >
Type deduction.
void build_sparsity_pattern(la::SparsityPattern &pattern, const Form< T, U > &a)
Build a sparsity pattern for a given form.
Definition utils.h:184
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
std::vector< std::int32_t > exterior_facet_indices(const Topology &topology, int facet_type_idx)
Compute the indices of all exterior facets that are owned by the caller.
Definition utils.cpp:58
Represents integral data, containing the kernel, and a list of entities to integrate over and the ind...
Definition Form.h:51