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/Topology.h>
24#include <dolfinx/mesh/cell_types.h>
25#include <dolfinx/mesh/utils.h>
26#include <functional>
27#include <memory>
28#include <optional>
29#include <set>
30#include <span>
31#include <stdexcept>
32#include <string>
33#include <tuple>
34#include <type_traits>
35#include <ufcx.h>
36#include <utility>
37#include <vector>
38
41namespace basix
42{
43template <std::floating_point T>
44class FiniteElement;
45}
46
47namespace dolfinx::common
48{
49class IndexMap;
50}
51
52namespace dolfinx::fem
53{
54template <dolfinx::scalar T, std::floating_point U>
55class Expression;
56
57namespace impl
58{
65template <int num_cells>
66std::array<std::int32_t, 2 * num_cells>
67get_cell_facet_pairs(std::int32_t f, std::span<const std::int32_t> cells,
69{
70 // Loop over cells sharing facet
71 assert(cells.size() == num_cells);
72 std::array<std::int32_t, 2 * num_cells> cell_local_facet_pairs;
73 for (int c = 0; c < num_cells; ++c)
74 {
75 // Get local index of facet with respect to the cell
76 std::int32_t cell = cells[c];
77 auto cell_facets = c_to_f.links(cell);
78 auto facet_it = std::find(cell_facets.begin(), cell_facets.end(), f);
79 assert(facet_it != cell_facets.end());
80 int local_f = std::distance(cell_facets.begin(), facet_it);
81 cell_local_facet_pairs[2 * c] = cell;
82 cell_local_facet_pairs[2 * c + 1] = local_f;
83 }
84
85 return cell_local_facet_pairs;
86}
87} // namespace impl
88
115std::vector<std::int32_t>
117 const mesh::Topology& topology,
118 std::span<const std::int32_t> entities);
119
127template <dolfinx::scalar T, std::floating_point U>
128std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>>
129extract_function_spaces(const std::vector<std::vector<const Form<T, U>*>>& a)
130{
131 std::vector<
132 std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>>
133 spaces(
134 a.size(),
135 std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>(
136 a.front().size()));
137 for (std::size_t i = 0; i < a.size(); ++i)
138 {
139 for (std::size_t j = 0; j < a[i].size(); ++j)
140 {
141 if (const Form<T, U>* form = a[i][j]; form)
142 spaces[i][j] = {form->function_spaces()[0], form->function_spaces()[1]};
143 }
144 }
145 return spaces;
146}
147
153template <dolfinx::scalar T, std::floating_point U>
155{
156 std::shared_ptr mesh = a.mesh();
157 assert(mesh);
158
159 // Get index maps and block sizes from the DOF maps. Note that in
160 // mixed-topology meshes, despite there being multiple DOF maps, the
161 // index maps and block sizes are the same.
162 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
163 *a.function_spaces().at(0)->dofmaps(0),
164 *a.function_spaces().at(1)->dofmaps(0)};
165
166 const std::array index_maps{dofmaps[0].get().index_map,
167 dofmaps[1].get().index_map};
168 const std::array bs
169 = {dofmaps[0].get().index_map_bs(), dofmaps[1].get().index_map_bs()};
170
171 la::SparsityPattern pattern(mesh->comm(), index_maps, bs);
172 build_sparsity_pattern(pattern, a);
173 return pattern;
174}
175
181template <dolfinx::scalar T, std::floating_point U>
183{
184 if (a.rank() != 2)
185 {
186 throw std::runtime_error(
187 "Cannot create sparsity pattern. Form is not a bilinear.");
188 }
189
190 std::shared_ptr mesh = a.mesh();
191 assert(mesh);
192 std::shared_ptr mesh0 = a.function_spaces().at(0)->mesh();
193 assert(mesh0);
194 std::shared_ptr mesh1 = a.function_spaces().at(1)->mesh();
195 assert(mesh1);
196
197 const std::set<IntegralType> types = a.integral_types();
198 if (types.find(IntegralType::interior_facet) != types.end()
199 or types.find(IntegralType::exterior_facet) != types.end())
200 {
201 // FIXME: cleanup these calls? Some of the happen internally again.
202 int tdim = mesh->topology()->dim();
203 mesh->topology_mutable()->create_entities(tdim - 1);
204 mesh->topology_mutable()->create_connectivity(tdim - 1, tdim);
205 }
206
207 common::Timer t0("Build sparsity");
208
209 auto extract_cells = [](std::span<const std::int32_t> facets)
210 {
211 assert(facets.size() % 2 == 0);
212 std::vector<std::int32_t> cells;
213 cells.reserve(facets.size() / 2);
214 for (std::size_t i = 0; i < facets.size(); i += 2)
215 cells.push_back(facets[i]);
216 return cells;
217 };
218
219 const int num_cell_types = mesh->topology()->cell_types().size();
220 for (int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
221 {
222 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
223 *a.function_spaces().at(0)->dofmaps(cell_type_idx),
224 *a.function_spaces().at(1)->dofmaps(cell_type_idx)};
225
226 // Create and build sparsity pattern
227 for (auto type : types)
228 {
229 std::vector<int> ids = a.integral_ids(type);
230 switch (type)
231 {
233 for (int id : ids)
234 {
235 sparsitybuild::cells(pattern,
236 {a.domain_arg(type, 0, id, cell_type_idx),
237 a.domain_arg(type, 1, id, cell_type_idx)},
238 {{dofmaps[0], dofmaps[1]}});
239 }
240 break;
242 for (int id : ids)
243 {
245 pattern,
246 {extract_cells(a.domain_arg(type, 0, id, 0)),
247 extract_cells(a.domain_arg(type, 1, id, 0))},
248 {{dofmaps[0], dofmaps[1]}});
249 }
250 break;
252 for (int id : ids)
253 {
254 sparsitybuild::cells(pattern,
255 {extract_cells(a.domain_arg(type, 0, id, 0)),
256 extract_cells(a.domain_arg(type, 1, id, 0))},
257 {{dofmaps[0], dofmaps[1]}});
258 }
259 break;
260 default:
261 throw std::runtime_error("Unsupported integral type");
262 }
263 }
264 }
265
266 t0.stop();
267}
268
270template <std::floating_point T>
272 const std::vector<int>& parent_map
273 = {})
274{
275 // Create subdofmaps and compute offset
276 std::vector<int> offsets(1, 0);
277 std::vector<dolfinx::fem::ElementDofLayout> sub_doflayout;
278 int bs = element.block_size();
279 for (int i = 0; i < element.num_sub_elements(); ++i)
280 {
281 // The ith sub-element. For mixed elements this is subelements()[i]. For
282 // blocked elements, the sub-element will always be the same, so we'll use
283 // sub_elements()[0]
284 std::shared_ptr<const fem::FiniteElement<T>> sub_e
285 = element.sub_elements()[bs > 1 ? 0 : i];
286
287 // In a mixed element DOFs are ordered element by element, so the offset to
288 // the next sub-element is sub_e->space_dimension(). Blocked elements use
289 // xxyyzz ordering, so the offset to the next sub-element is 1
290
291 std::vector<int> parent_map_sub(sub_e->space_dimension(), offsets.back());
292 for (std::size_t j = 0; j < parent_map_sub.size(); ++j)
293 parent_map_sub[j] += bs * j;
294 offsets.push_back(offsets.back() + (bs > 1 ? 1 : sub_e->space_dimension()));
295 sub_doflayout.push_back(
296 dolfinx::fem::create_element_dof_layout(*sub_e, parent_map_sub));
297 }
298
299 return ElementDofLayout(bs, element.entity_dofs(),
300 element.entity_closure_dofs(), parent_map,
301 sub_doflayout);
302}
303
313 MPI_Comm comm, const ElementDofLayout& layout, mesh::Topology& topology,
314 std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv,
315 std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
316 reorder_fn);
317
328std::vector<DofMap> create_dofmaps(
329 MPI_Comm comm, const std::vector<ElementDofLayout>& layouts,
330 mesh::Topology& topology,
331 std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv,
332 std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
333 reorder_fn);
334
338std::vector<std::string> get_coefficient_names(const ufcx_form& ufcx_form);
339
343std::vector<std::string> get_constant_names(const ufcx_form& ufcx_form);
344
363template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
365 const std::vector<std::reference_wrapper<const ufcx_form>>& ufcx_forms,
366 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
367 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients,
368 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
369 const std::map<
371 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
372 subdomains,
373 const std::map<std::shared_ptr<const mesh::Mesh<U>>,
374 std::span<const std::int32_t>>& entity_maps,
375 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
376{
377 for (const ufcx_form& ufcx_form : ufcx_forms)
378 {
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())
382 {
383 throw std::runtime_error("Mismatch between number of expected and "
384 "provided Form coefficients.");
385 }
386 if (ufcx_form.num_constants != (int)constants.size())
387 {
388 throw std::runtime_error(
389 "Mismatch between number of expected and provided Form constants.");
390 }
391 }
392
393 // Check argument function spaces
394 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
395 {
396 for (std::size_t i = 0; i < spaces.size(); ++i)
397 {
398 assert(spaces[i]->elements(form_idx));
399 if (auto element_hash
400 = ufcx_forms[form_idx].get().finite_element_hashes[i];
401 element_hash != 0
402 and element_hash
403 != spaces[i]->elements(form_idx)->basix_element().hash())
404 {
405 throw std::runtime_error(
406 "Cannot create form. Elements are different to "
407 "those used to compile the form.");
408 }
409 }
410 }
411
412 // Extract mesh from FunctionSpace, and check they are the same
413 if (!mesh and !spaces.empty())
414 mesh = spaces.front()->mesh();
415 if (!mesh)
416 throw std::runtime_error("No mesh could be associated with the Form.");
417 for (auto& V : spaces)
418 {
419 if (mesh != V->mesh() and entity_maps.find(V->mesh()) == entity_maps.end())
420 {
421 throw std::runtime_error(
422 "Incompatible mesh. entity_maps must be provided.");
423 }
424 }
425
426 auto topology = mesh->topology();
427 assert(topology);
428 const int tdim = topology->dim();
429
430 // NOTE: This assumes all forms in mixed-topology meshes have the same
431 // integral offsets. Since the UFL forms for each type of cell should be
432 // the same, I think this assumption is OK.
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];
437
438 // Create facets, if required
439 if (num_integrals_type[exterior_facet] > 0
440 or num_integrals_type[interior_facet] > 0)
441 {
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);
445 }
446
447 // Get list of integral IDs, and load tabulate tensor into memory for
448 // each
449 using kern_t = std::function<void(T*, const T*, const T*, const U*,
450 const int*, const std::uint8_t*, void*)>;
451 std::map<std::tuple<IntegralType, int, int>, integral_data<T, U>> integrals;
452
453 auto check_geometry_hash
454 = [&geo = mesh->geometry()](const ufcx_integral& integral,
455 std::size_t cell_idx)
456 {
457 if (integral.coordinate_element_hash != geo.cmaps().at(cell_idx).hash())
458 {
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()));
463 }
464 };
465
466 // Attach cell kernels
467 bool needs_facet_permutations = false;
468 {
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]);
473 auto sd = subdomains.find(IntegralType::cell);
474 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
475 {
476 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
477 for (int i = 0; i < num_integrals_type[cell]; ++i)
478 {
479 const int id = ids[i];
480 ufcx_integral* integral
481 = ufcx_form.form_integrals[integral_offsets[cell] + i];
482 assert(integral);
483 check_geometry_hash(*integral, form_idx);
484
485 // Build list of active coefficients
486 std::vector<int> active_coeffs;
487 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
488 {
489 if (integral->enabled_coefficients[j])
490 active_coeffs.push_back(j);
491 }
492
493 kern_t k = nullptr;
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>>)
498 {
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);
503 }
504#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
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>>)
509 {
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);
514 }
515#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
516
517 if (!k)
518 {
519 throw std::runtime_error(
520 "UFCx kernel function is NULL. Check requested types.");
521 }
522
523 // Build list of entities to assemble over
524 if (id == -1)
525 {
526 // Default kernel, operates on all (owned) cells
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);
531 integrals.insert({{IntegralType::cell, id, form_idx},
532 {k, default_cells, active_coeffs}});
533 }
534 else if (sd != subdomains.end())
535 {
536 // NOTE: This requires that pairs are sorted
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)
540 {
541 integrals.insert({{IntegralType::cell, id, form_idx},
542 {k,
543 std::vector<std::int32_t>(it->second.begin(),
544 it->second.end()),
545 active_coeffs}});
546 }
547 }
548
549 if (integral->needs_facet_permutations)
550 needs_facet_permutations = true;
551 }
552 }
553 }
554
555 // Attach exterior facet kernels
556 std::vector<std::int32_t> default_facets_ext;
557 {
558 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
559 + integral_offsets[exterior_facet],
560 num_integrals_type[exterior_facet]);
561 auto sd = subdomains.find(IntegralType::exterior_facet);
562 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
563 {
564 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
565 for (int i = 0; i < num_integrals_type[exterior_facet]; ++i)
566 {
567 const int id = ids[i];
568 ufcx_integral* integral
569 = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i];
570 assert(integral);
571 check_geometry_hash(*integral, form_idx);
572
573 std::vector<int> active_coeffs;
574 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
575 {
576 if (integral->enabled_coefficients[j])
577 active_coeffs.push_back(j);
578 }
579
580 kern_t k = nullptr;
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>>)
585 {
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);
590 }
591#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
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>>)
596 {
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);
601 }
602#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
603 assert(k);
604
605 // Build list of entities to assembler over
606 const std::vector bfacets = mesh::exterior_facet_indices(*topology);
607 auto f_to_c = topology->connectivity(tdim - 1, tdim);
608 assert(f_to_c);
609 auto c_to_f = topology->connectivity(tdim, tdim - 1);
610 assert(c_to_f);
611 if (id == -1)
612 {
613 // Default kernel, operates on all (owned) exterior facets
614 default_facets_ext.reserve(2 * bfacets.size());
615 for (std::int32_t f : bfacets)
616 {
617 // There will only be one pair for an exterior facet integral
618 auto pair
619 = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f);
620 default_facets_ext.insert(default_facets_ext.end(), pair.begin(),
621 pair.end());
622 }
623 integrals.insert({{IntegralType::exterior_facet, id, form_idx},
624 {k, default_facets_ext, active_coeffs}});
625 }
626 else if (sd != subdomains.end())
627 {
628 // NOTE: This requires that pairs are sorted
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)
632 {
633 integrals.insert({{IntegralType::exterior_facet, id, form_idx},
634 {k,
635 std::vector<std::int32_t>(it->second.begin(),
636 it->second.end()),
637 active_coeffs}});
638 }
639 }
640
641 if (integral->needs_facet_permutations)
642 needs_facet_permutations = true;
643 }
644 }
645 }
646
647 // Attach interior facet kernels
648 std::vector<std::int32_t> default_facets_int;
649 {
650 std::span<const int> ids(ufcx_forms[0].get().form_integral_ids
651 + integral_offsets[interior_facet],
652 num_integrals_type[interior_facet]);
653 auto sd = subdomains.find(IntegralType::interior_facet);
654 for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx)
655 {
656 const ufcx_form& ufcx_form = ufcx_forms[form_idx];
657
658 // Create indicator for interprocess facets
659 std::vector<std::int8_t> interprocess_marker;
660 if (num_integrals_type[interior_facet] > 0)
661 {
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; });
671 }
672
673 for (int i = 0; i < num_integrals_type[interior_facet]; ++i)
674 {
675 const int id = ids[i];
676 ufcx_integral* integral
677 = ufcx_form.form_integrals[integral_offsets[interior_facet] + i];
678 assert(integral);
679 check_geometry_hash(*integral, form_idx);
680
681 std::vector<int> active_coeffs;
682 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
683 {
684 if (integral->enabled_coefficients[j])
685 active_coeffs.push_back(j);
686 }
687
688 kern_t k = nullptr;
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>>)
693 {
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);
698 }
699#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
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>>)
704 {
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);
709 }
710#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
711 assert(k);
712
713 // Build list of entities to assembler over
714 auto f_to_c = topology->connectivity(tdim - 1, tdim);
715 assert(f_to_c);
716 auto c_to_f = topology->connectivity(tdim, tdim - 1);
717 assert(c_to_f);
718 if (id == -1)
719 {
720 // Default kernel, operates on all (owned) interior facets
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)
725 {
726 if (f_to_c->num_links(f) == 2)
727 {
728 auto pairs
729 = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f);
730 default_facets_int.insert(default_facets_int.end(), pairs.begin(),
731 pairs.end());
732 }
733 else if (interprocess_marker[f])
734 {
735 throw std::runtime_error(
736 "Cannot compute interior facet integral over interprocess "
737 "facet. Please use ghost mode shared facet when creating the "
738 "mesh");
739 }
740 }
741 integrals.insert({{IntegralType::interior_facet, id, form_idx},
742 {k, default_facets_int, active_coeffs}});
743 }
744 else if (sd != subdomains.end())
745 {
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)
749 {
750 integrals.insert({{IntegralType::interior_facet, id, form_idx},
751 {k,
752 std::vector<std::int32_t>(it->second.begin(),
753 it->second.end()),
754 active_coeffs}});
755 }
756 }
757
758 if (integral->needs_facet_permutations)
759 needs_facet_permutations = true;
760 }
761 }
762 }
763
764 return Form<T, U>(spaces, std::move(integrals), mesh, coefficients, constants,
765 needs_facet_permutations, entity_maps);
766}
767
782template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
784 const ufcx_form& ufcx_form,
785 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
786 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
787 coefficients,
788 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
789 const std::map<
791 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
792 subdomains,
793 const std::map<std::shared_ptr<const mesh::Mesh<U>>,
794 std::span<const std::int32_t>>& entity_maps,
795 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
796{
797 // Place coefficients in appropriate order
798 std::vector<std::shared_ptr<const Function<T, U>>> coeff_map;
799 for (const std::string& name : get_coefficient_names(ufcx_form))
800 {
801 if (auto it = coefficients.find(name); it != coefficients.end())
802 coeff_map.push_back(it->second);
803 else
804 {
805 throw std::runtime_error("Form coefficient \"" + name
806 + "\" not provided.");
807 }
808 }
809
810 // Place constants in appropriate order
811 std::vector<std::shared_ptr<const Constant<T>>> const_map;
812 for (const std::string& name : get_constant_names(ufcx_form))
813 {
814 if (auto it = constants.find(name); it != constants.end())
815 const_map.push_back(it->second);
816 else
817 throw std::runtime_error("Form constant \"" + name + "\" not provided.");
818 }
819
820 return create_form_factory({ufcx_form}, spaces, coeff_map, const_map,
821 subdomains, entity_maps, mesh);
822}
823
842template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
844 ufcx_form* (*fptr)(),
845 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
846 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
847 coefficients,
848 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
849 const std::map<
851 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
852 subdomains,
853 const std::map<std::shared_ptr<const mesh::Mesh<U>>,
854 std::span<const std::int32_t>>& entity_maps,
855 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
856{
857 ufcx_form* form = fptr();
858 Form<T, U> L = create_form<T, U>(*form, spaces, coefficients, constants,
859 subdomains, entity_maps, mesh);
860 std::free(form);
861 return L;
862}
863
865template <std::floating_point T>
867 std::shared_ptr<mesh::Mesh<T>> mesh,
868 std::shared_ptr<const fem::FiniteElement<T>> e,
869 std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
870 reorder_fn
871 = nullptr)
872{
873 // TODO: check cell type of e (need to add method to fem::FiniteElement)
874 assert(e);
875 assert(mesh);
876 assert(mesh->topology());
877 if (e->cell_type() != mesh->topology()->cell_type())
878 throw std::runtime_error("Cell type of element and mesh must match.");
879
880 // Create element dof layout
882
883 // Create a dofmap
884 std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv
885 = e->needs_dof_permutations() ? e->dof_permutation_fn(true, true)
886 : nullptr;
887 auto dofmap = std::make_shared<const DofMap>(create_dofmap(
888 mesh->comm(), layout, *mesh->topology(), permute_inv, reorder_fn));
889
890 return FunctionSpace(mesh, e, dofmap);
891}
892
894template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
896 const ufcx_expression& e,
897 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients,
898 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
899 std::shared_ptr<const FunctionSpace<U>> argument_space = nullptr)
900{
901 if (!coefficients.empty())
902 {
903 assert(coefficients.front());
904 assert(coefficients.front()->function_space());
905 std::shared_ptr<const mesh::Mesh<U>> mesh
906 = coefficients.front()->function_space()->mesh();
907 if (mesh->geometry().cmap().hash() != e.coordinate_element_hash)
908 {
909 throw std::runtime_error(
910 "Expression and mesh geometric maps do not match.");
911 }
912 }
913
914 if (e.rank > 0 and !argument_space)
915 {
916 throw std::runtime_error("Expression has Argument but no Argument "
917 "function space was provided.");
918 }
919
920 std::vector<U> X(e.points, e.points + e.num_points * e.entity_dimension);
921 std::array<std::size_t, 2> Xshape
922 = {static_cast<std::size_t>(e.num_points),
923 static_cast<std::size_t>(e.entity_dimension)};
924 std::vector<std::size_t> value_shape(e.value_shape,
925 e.value_shape + e.num_components);
926 std::function<void(T*, const T*, const T*, const scalar_value_t<T>*,
927 const int*, const std::uint8_t*, void*)>
928 tabulate_tensor = nullptr;
929 if constexpr (std::is_same_v<T, float>)
930 tabulate_tensor = e.tabulate_tensor_float32;
931#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
932 else if constexpr (std::is_same_v<T, std::complex<float>>)
933 {
934 tabulate_tensor = reinterpret_cast<void (*)(
935 T*, const T*, const T*, const scalar_value_t<T>*, const int*,
936 const unsigned char*, void*)>(e.tabulate_tensor_complex64);
937 }
938#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
939 else if constexpr (std::is_same_v<T, double>)
940 tabulate_tensor = e.tabulate_tensor_float64;
941#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
942 else if constexpr (std::is_same_v<T, std::complex<double>>)
943 {
944 tabulate_tensor = reinterpret_cast<void (*)(
945 T*, const T*, const T*, const scalar_value_t<T>*, const int*,
946 const unsigned char*, void*)>(e.tabulate_tensor_complex128);
947 }
948#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
949 else
950 throw std::runtime_error("Type not supported.");
951
952 assert(tabulate_tensor);
953 return Expression(coefficients, constants, std::span<const U>(X), Xshape,
954 tabulate_tensor, value_shape, argument_space);
955}
956
959template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
961 const ufcx_expression& e,
962 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
963 coefficients,
964 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
965 std::shared_ptr<const FunctionSpace<U>> argument_space = nullptr)
966{
967 // Place coefficients in appropriate order
968 std::vector<std::shared_ptr<const Function<T, U>>> coeff_map;
969 std::vector<std::string> coefficient_names;
970 for (int i = 0; i < e.num_coefficients; ++i)
971 coefficient_names.push_back(e.coefficient_names[i]);
972
973 for (const std::string& name : coefficient_names)
974 {
975 if (auto it = coefficients.find(name); it != coefficients.end())
976 coeff_map.push_back(it->second);
977 else
978 {
979 throw std::runtime_error("Expression coefficient \"" + name
980 + "\" not provided.");
981 }
982 }
983
984 // Place constants in appropriate order
985 std::vector<std::shared_ptr<const Constant<T>>> const_map;
986 std::vector<std::string> constant_names;
987 for (int i = 0; i < e.num_constants; ++i)
988 constant_names.push_back(e.constant_names[i]);
989
990 for (const std::string& name : constant_names)
991 {
992 if (auto it = constants.find(name); it != constants.end())
993 const_map.push_back(it->second);
994 else
995 {
996 throw std::runtime_error("Expression constant \"" + name
997 + "\" not provided.");
998 }
999 }
1000
1001 return create_expression(e, coeff_map, const_map, argument_space);
1002}
1003
1004} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
Definition CoordinateElement.h:26
Definition IndexMap.h:94
Timer for measuring and logging elapsed time durations.
Definition Timer.h:41
std::chrono::duration< double, Period > stop()
Stop timer and return elapsed time.
Definition Timer.h:91
Constant value which can be attached to a Form.
Definition Constant.h:23
Degree-of-freedom map.
Definition DofMap.h:73
Definition ElementDofLayout.h:30
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:175
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
Definition Function.h:47
Definition AdjacencyList.h:27
std::span< T > links(std::size_t node)
Definition AdjacencyList.h:111
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:46
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:67
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:124
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:866
std::vector< DofMap > create_dofmaps(MPI_Comm comm, const std::vector< ElementDofLayout > &layouts, mesh::Topology &topology, std::function< void(std::span< std::int32_t >, std::uint32_t)> permute_inv, 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:68
DofMap create_dofmap(MPI_Comm comm, const ElementDofLayout &layout, mesh::Topology &topology, std::function< void(std::span< std::int32_t >, std::uint32_t)> permute_inv, 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< std::string > get_coefficient_names(const ufcx_form &ufcx_form)
Definition utils.cpp:117
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:129
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:895
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:132
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::map< std::shared_ptr< const mesh::Mesh< U > >, std::span< const std::int32_t > > &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:364
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:271
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::map< std::shared_ptr< const mesh::Mesh< U > >, std::span< const std::int32_t > > &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:783
IntegralType
Type of integral.
Definition Form.h:36
@ interior_facet
Interior facet.
Definition Form.h:39
@ cell
Cell.
Definition Form.h:37
@ exterior_facet
Exterior facet.
Definition Form.h:38
la::SparsityPattern create_sparsity_pattern(const Form< T, U > &a)
Create a sparsity pattern for a given form.
Definition utils.h:154
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:182
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:110