DOLFINx 0.9.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 <set>
29#include <span>
30#include <stdexcept>
31#include <string>
32#include <type_traits>
33#include <ufcx.h>
34#include <utility>
35#include <vector>
36
39namespace basix
40{
41template <std::floating_point T>
43}
44
45namespace dolfinx::common
46{
47class IndexMap;
48}
49
50namespace dolfinx::fem
51{
52
53namespace impl
54{
61template <int num_cells>
62std::array<std::int32_t, 2 * num_cells>
63get_cell_facet_pairs(std::int32_t f, std::span<const std::int32_t> cells,
65{
66 // Loop over cells sharing facet
67 assert(cells.size() == num_cells);
68 std::array<std::int32_t, 2 * num_cells> cell_local_facet_pairs;
69 for (int c = 0; c < num_cells; ++c)
70 {
71 // Get local index of facet with respect to the cell
72 std::int32_t cell = cells[c];
73 auto cell_facets = c_to_f.links(cell);
74 auto facet_it = std::find(cell_facets.begin(), cell_facets.end(), f);
75 assert(facet_it != cell_facets.end());
76 int local_f = std::distance(cell_facets.begin(), facet_it);
77 cell_local_facet_pairs[2 * c] = cell;
78 cell_local_facet_pairs[2 * c + 1] = local_f;
79 }
80
81 return cell_local_facet_pairs;
82}
83
84} // namespace impl
85
108std::vector<std::int32_t>
110 const mesh::Topology& topology,
111 std::span<const std::int32_t> entities, int dim);
112
120template <dolfinx::scalar T, std::floating_point U>
121std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>>
122extract_function_spaces(const std::vector<std::vector<const Form<T, U>*>>& a)
123{
124 std::vector<
125 std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>>
126 spaces(
127 a.size(),
128 std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 2>>(
129 a.front().size()));
130 for (std::size_t i = 0; i < a.size(); ++i)
131 {
132 for (std::size_t j = 0; j < a[i].size(); ++j)
133 {
134 if (const Form<T, U>* form = a[i][j]; form)
135 spaces[i][j] = {form->function_spaces()[0], form->function_spaces()[1]};
136 }
137 }
138 return spaces;
139}
140
146template <dolfinx::scalar T, std::floating_point U>
148{
149 if (a.rank() != 2)
150 {
151 throw std::runtime_error(
152 "Cannot create sparsity pattern. Form is not a bilinear.");
153 }
154
155 // Get dof maps and mesh
156 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
157 *a.function_spaces().at(0)->dofmap(),
158 *a.function_spaces().at(1)->dofmap()};
159 std::shared_ptr mesh = a.mesh();
160 assert(mesh);
161
162 std::shared_ptr mesh0 = a.function_spaces().at(0)->mesh();
163 assert(mesh0);
164 std::shared_ptr mesh1 = a.function_spaces().at(1)->mesh();
165 assert(mesh1);
166
167 const std::set<IntegralType> types = a.integral_types();
168 if (types.find(IntegralType::interior_facet) != types.end()
169 or types.find(IntegralType::exterior_facet) != types.end())
170 {
171 // FIXME: cleanup these calls? Some of the happen internally again.
172 int tdim = mesh->topology()->dim();
173 mesh->topology_mutable()->create_entities(tdim - 1);
174 mesh->topology_mutable()->create_connectivity(tdim - 1, tdim);
175 }
176
177 common::Timer t0("Build sparsity");
178
179 // Get common::IndexMaps for each dimension
180 const std::array index_maps{dofmaps[0].get().index_map,
181 dofmaps[1].get().index_map};
182 const std::array bs
183 = {dofmaps[0].get().index_map_bs(), dofmaps[1].get().index_map_bs()};
184
185 auto extract_cells = [](std::span<const std::int32_t> facets)
186 {
187 assert(facets.size() % 2 == 0);
188 std::vector<std::int32_t> cells;
189 cells.reserve(facets.size() / 2);
190 for (std::size_t i = 0; i < facets.size(); i += 2)
191 cells.push_back(facets[i]);
192 return cells;
193 };
194
195 // Create and build sparsity pattern
196 la::SparsityPattern pattern(mesh->comm(), index_maps, bs);
197 for (auto type : types)
198 {
199 std::vector<int> ids = a.integral_ids(type);
200 switch (type)
201 {
203 for (int id : ids)
204 {
206 pattern, {a.domain(type, id, *mesh0), a.domain(type, id, *mesh1)},
207 {{dofmaps[0], dofmaps[1]}});
208 }
209 break;
211 for (int id : ids)
212 {
214 pattern,
215 {extract_cells(a.domain(type, id, *mesh0)),
216 extract_cells(a.domain(type, id, *mesh1))},
217 {{dofmaps[0], dofmaps[1]}});
218 }
219 break;
221 for (int id : ids)
222 {
223 sparsitybuild::cells(pattern,
224 {extract_cells(a.domain(type, id, *mesh0)),
225 extract_cells(a.domain(type, id, *mesh1))},
226 {{dofmaps[0], dofmaps[1]}});
227 }
228 break;
229 default:
230 throw std::runtime_error("Unsupported integral type");
231 }
232 }
233
234 t0.stop();
235
236 return pattern;
237}
238
240template <std::floating_point T>
242 const std::vector<int>& parent_map
243 = {})
244{
245 // Create subdofmaps and compute offset
246 std::vector<int> offsets(1, 0);
247 std::vector<dolfinx::fem::ElementDofLayout> sub_doflayout;
248 int bs = element.block_size();
249 for (int i = 0; i < element.num_sub_elements(); ++i)
250 {
251 // The ith sub-element. For mixed elements this is subelements()[i]. For
252 // blocked elements, the sub-element will always be the same, so we'll use
253 // sub_elements()[0]
254 std::shared_ptr<const fem::FiniteElement<T>> sub_e
255 = element.sub_elements()[bs > 1 ? 0 : i];
256
257 // In a mixed element DOFs are ordered element by element, so the offset to
258 // the next sub-element is sub_e->space_dimension(). Blocked elements use
259 // xxyyzz ordering, so the offset to the next sub-element is 1
260
261 std::vector<int> parent_map_sub(sub_e->space_dimension(), offsets.back());
262 for (std::size_t j = 0; j < parent_map_sub.size(); ++j)
263 parent_map_sub[j] += bs * j;
264 offsets.push_back(offsets.back() + (bs > 1 ? 1 : sub_e->space_dimension()));
265 sub_doflayout.push_back(
266 dolfinx::fem::create_element_dof_layout(*sub_e, parent_map_sub));
267 }
268
269 return ElementDofLayout(bs, element.entity_dofs(),
270 element.entity_closure_dofs(), parent_map,
271 sub_doflayout);
272}
273
282DofMap create_dofmap(
283 MPI_Comm comm, const ElementDofLayout& layout, mesh::Topology& topology,
284 std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv,
285 std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
286 reorder_fn);
287
298std::vector<DofMap> create_dofmaps(
299 MPI_Comm comm, const std::vector<ElementDofLayout>& layouts,
300 mesh::Topology& topology,
301 std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv,
302 std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
303 reorder_fn);
304
308std::vector<std::string> get_coefficient_names(const ufcx_form& ufcx_form);
309
313std::vector<std::string> get_constant_names(const ufcx_form& ufcx_form);
314
332template <dolfinx::scalar T, std::floating_point U = scalar_value_type_t<T>>
334 const ufcx_form& ufcx_form,
335 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
336 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients,
337 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
338 const std::map<
340 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
341 subdomains,
342 const std::map<std::shared_ptr<const mesh::Mesh<U>>,
343 std::span<const std::int32_t>>& entity_maps,
344 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
345{
346 if (ufcx_form.rank != (int)spaces.size())
347 throw std::runtime_error("Wrong number of argument spaces for Form.");
348 if (ufcx_form.num_coefficients != (int)coefficients.size())
349 {
350 throw std::runtime_error(
351 "Mismatch between number of expected and provided Form coefficients.");
352 }
353 if (ufcx_form.num_constants != (int)constants.size())
354 {
355 throw std::runtime_error(
356 "Mismatch between number of expected and provided Form constants.");
357 }
358
359 // Check argument function spaces
360 for (std::size_t i = 0; i < spaces.size(); ++i)
361 {
362 assert(spaces[i]->element());
363 if (auto element_hash = ufcx_form.finite_element_hashes[i];
364 element_hash != 0
365 and element_hash != spaces[i]->element()->basix_element().hash())
366 {
367 throw std::runtime_error("Cannot create form. Elements are different to "
368 "those used to compile the form.");
369 }
370 }
371
372 // Extract mesh from FunctionSpace, and check they are the same
373 if (!mesh and !spaces.empty())
374 mesh = spaces[0]->mesh();
375 for (auto& V : spaces)
376 {
377 if (mesh != V->mesh() and entity_maps.find(V->mesh()) == entity_maps.end())
378 throw std::runtime_error(
379 "Incompatible mesh. entity_maps must be provided.");
380 }
381 if (!mesh)
382 throw std::runtime_error("No mesh could be associated with the Form.");
383
384 auto topology = mesh->topology();
385 assert(topology);
386 const int tdim = topology->dim();
387
388 const int* integral_offsets = ufcx_form.form_integral_offsets;
389 std::vector<int> num_integrals_type(3);
390 for (int i = 0; i < 3; ++i)
391 num_integrals_type[i] = integral_offsets[i + 1] - integral_offsets[i];
392
393 // Create facets, if required
394 if (num_integrals_type[exterior_facet] > 0
395 or num_integrals_type[interior_facet] > 0)
396 {
397 mesh->topology_mutable()->create_entities(tdim - 1);
398 mesh->topology_mutable()->create_connectivity(tdim - 1, tdim);
399 mesh->topology_mutable()->create_connectivity(tdim, tdim - 1);
400 }
401
402 // Get list of integral IDs, and load tabulate tensor into memory for
403 // each
404 using kern_t = std::function<void(T*, const T*, const T*, const U*,
405 const int*, const std::uint8_t*)>;
406 std::map<IntegralType, std::vector<integral_data<T, U>>> integrals;
407
408 // Attach cell kernels
409 bool needs_facet_permutations = false;
410 std::vector<std::int32_t> default_cells;
411 {
412 std::span<const int> ids(ufcx_form.form_integral_ids
413 + integral_offsets[cell],
414 num_integrals_type[cell]);
415 auto itg = integrals.insert({IntegralType::cell, {}});
416 auto sd = subdomains.find(IntegralType::cell);
417 for (int i = 0; i < num_integrals_type[cell]; ++i)
418 {
419 const int id = ids[i];
420 ufcx_integral* integral
421 = ufcx_form.form_integrals[integral_offsets[cell] + i];
422 assert(integral);
423
424 // Build list of active coefficients
425 std::vector<int> active_coeffs;
426 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
427 {
428 if (integral->enabled_coefficients[j])
429 active_coeffs.push_back(j);
430 }
431
432 kern_t k = nullptr;
433 if constexpr (std::is_same_v<T, float>)
434 k = integral->tabulate_tensor_float32;
435#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
436 else if constexpr (std::is_same_v<T, std::complex<float>>)
437 {
438 k = reinterpret_cast<void (*)(
439 T*, const T*, const T*,
440 const typename scalar_value_type<T>::value_type*, const int*,
441 const unsigned char*)>(integral->tabulate_tensor_complex64);
442 }
443#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
444 else if constexpr (std::is_same_v<T, double>)
445 k = integral->tabulate_tensor_float64;
446#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
447 else if constexpr (std::is_same_v<T, std::complex<double>>)
448 {
449 k = reinterpret_cast<void (*)(
450 T*, const T*, const T*,
451 const typename scalar_value_type<T>::value_type*, const int*,
452 const unsigned char*)>(integral->tabulate_tensor_complex128);
453 }
454#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
455
456 if (!k)
457 {
458 throw std::runtime_error(
459 "UFCx kernel function is NULL. Check requested types.");
460 }
461
462 // Build list of entities to assemble over
463 if (id == -1)
464 {
465 // Default kernel, operates on all (owned) cells
466 assert(topology->index_map(tdim));
467 default_cells.resize(topology->index_map(tdim)->size_local(), 0);
468 std::iota(default_cells.begin(), default_cells.end(), 0);
469 itg.first->second.emplace_back(id, k, default_cells, active_coeffs);
470 }
471 else if (sd != subdomains.end())
472 {
473 // NOTE: This requires that pairs are sorted
474 auto it
475 = std::ranges::lower_bound(sd->second, id, std::less<>{},
476 [](const auto& a) { return a.first; });
477 if (it != sd->second.end() and it->first == id)
478 itg.first->second.emplace_back(id, k, it->second, active_coeffs);
479 }
480
481 if (integral->needs_facet_permutations)
482 needs_facet_permutations = true;
483 }
484 }
485
486 // Attach exterior facet kernels
487 std::vector<std::int32_t> default_facets_ext;
488 {
489 std::span<const int> ids(ufcx_form.form_integral_ids
490 + integral_offsets[exterior_facet],
491 num_integrals_type[exterior_facet]);
492 auto itg = integrals.insert({IntegralType::exterior_facet, {}});
493 auto sd = subdomains.find(IntegralType::exterior_facet);
494 for (int i = 0; i < num_integrals_type[exterior_facet]; ++i)
495 {
496 const int id = ids[i];
497 ufcx_integral* integral
498 = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i];
499 assert(integral);
500 std::vector<int> active_coeffs;
501 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
502 {
503 if (integral->enabled_coefficients[j])
504 active_coeffs.push_back(j);
505 }
506
507 kern_t k = nullptr;
508 if constexpr (std::is_same_v<T, float>)
509 k = integral->tabulate_tensor_float32;
510#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
511 else if constexpr (std::is_same_v<T, std::complex<float>>)
512 {
513 k = reinterpret_cast<void (*)(
514 T*, const T*, const T*,
515 const typename scalar_value_type<T>::value_type*, const int*,
516 const unsigned char*)>(integral->tabulate_tensor_complex64);
517 }
518#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
519 else if constexpr (std::is_same_v<T, double>)
520 k = integral->tabulate_tensor_float64;
521#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
522 else if constexpr (std::is_same_v<T, std::complex<double>>)
523 {
524 k = reinterpret_cast<void (*)(
525 T*, const T*, const T*,
526 const typename scalar_value_type<T>::value_type*, const int*,
527 const unsigned char*)>(integral->tabulate_tensor_complex128);
528 }
529#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
530 assert(k);
531
532 // Build list of entities to assembler over
533 const std::vector bfacets = mesh::exterior_facet_indices(*topology);
534 auto f_to_c = topology->connectivity(tdim - 1, tdim);
535 assert(f_to_c);
536 auto c_to_f = topology->connectivity(tdim, tdim - 1);
537 assert(c_to_f);
538 if (id == -1)
539 {
540 // Default kernel, operates on all (owned) exterior facets
541 default_facets_ext.reserve(2 * bfacets.size());
542 for (std::int32_t f : bfacets)
543 {
544 // There will only be one pair for an exterior facet integral
545 auto pair
546 = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f);
547 default_facets_ext.insert(default_facets_ext.end(), pair.begin(),
548 pair.end());
549 }
550 itg.first->second.emplace_back(id, k, default_facets_ext,
551 active_coeffs);
552 }
553 else if (sd != subdomains.end())
554 {
555 // NOTE: This requires that pairs are sorted
556 auto it
557 = std::ranges::lower_bound(sd->second, id, std::less<>{},
558 [](const auto& a) { return a.first; });
559 if (it != sd->second.end() and it->first == id)
560 itg.first->second.emplace_back(id, k, it->second, active_coeffs);
561 }
562
563 if (integral->needs_facet_permutations)
564 needs_facet_permutations = true;
565 }
566 }
567
568 // Attach interior facet kernels
569 std::vector<std::int32_t> default_facets_int;
570 {
571 std::span<const int> ids(ufcx_form.form_integral_ids
572 + integral_offsets[interior_facet],
573 num_integrals_type[interior_facet]);
574 auto itg = integrals.insert({IntegralType::interior_facet, {}});
575 auto sd = subdomains.find(IntegralType::interior_facet);
576
577 // Create indicator for interprocess facets
578 std::vector<std::int8_t> interprocess_marker;
579 if (num_integrals_type[interior_facet] > 0)
580 {
581 assert(topology->index_map(tdim - 1));
582 const std::vector<std::int32_t>& interprocess_facets
583 = topology->interprocess_facets();
584 std::int32_t num_facets = topology->index_map(tdim - 1)->size_local()
585 + topology->index_map(tdim - 1)->num_ghosts();
586 interprocess_marker.resize(num_facets, 0);
587 std::ranges::for_each(interprocess_facets, [&interprocess_marker](auto f)
588 { interprocess_marker[f] = 1; });
589 }
590
591 for (int i = 0; i < num_integrals_type[interior_facet]; ++i)
592 {
593 const int id = ids[i];
594 ufcx_integral* integral
595 = ufcx_form.form_integrals[integral_offsets[interior_facet] + i];
596 assert(integral);
597 std::vector<int> active_coeffs;
598 for (int j = 0; j < ufcx_form.num_coefficients; ++j)
599 {
600 if (integral->enabled_coefficients[j])
601 active_coeffs.push_back(j);
602 }
603
604 kern_t k = nullptr;
605 if constexpr (std::is_same_v<T, float>)
606 k = integral->tabulate_tensor_float32;
607#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
608 else if constexpr (std::is_same_v<T, std::complex<float>>)
609 {
610 k = reinterpret_cast<void (*)(
611 T*, const T*, const T*,
612 const typename scalar_value_type<T>::value_type*, const int*,
613 const unsigned char*)>(integral->tabulate_tensor_complex64);
614 }
615#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
616 else if constexpr (std::is_same_v<T, double>)
617 k = integral->tabulate_tensor_float64;
618#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
619 else if constexpr (std::is_same_v<T, std::complex<double>>)
620 {
621 k = reinterpret_cast<void (*)(
622 T*, const T*, const T*,
623 const typename scalar_value_type<T>::value_type*, const int*,
624 const unsigned char*)>(integral->tabulate_tensor_complex128);
625 }
626#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
627 assert(k);
628
629 // Build list of entities to assembler over
630 auto f_to_c = topology->connectivity(tdim - 1, tdim);
631 assert(f_to_c);
632 auto c_to_f = topology->connectivity(tdim, tdim - 1);
633 assert(c_to_f);
634 if (id == -1)
635 {
636 // Default kernel, operates on all (owned) interior facets
637 assert(topology->index_map(tdim - 1));
638 std::int32_t num_facets = topology->index_map(tdim - 1)->size_local();
639 default_facets_int.reserve(4 * num_facets);
640 for (std::int32_t f = 0; f < num_facets; ++f)
641 {
642 if (f_to_c->num_links(f) == 2)
643 {
644 auto pairs
645 = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f);
646 default_facets_int.insert(default_facets_int.end(), pairs.begin(),
647 pairs.end());
648 }
649 else if (interprocess_marker[f])
650 {
651 throw std::runtime_error(
652 "Cannot compute interior facet integral over interprocess "
653 "facet. Please use ghost mode shared facet when creating the "
654 "mesh");
655 }
656 }
657 itg.first->second.emplace_back(id, k, default_facets_int,
658 active_coeffs);
659 }
660 else if (sd != subdomains.end())
661 {
662 auto it
663 = std::ranges::lower_bound(sd->second, id, std::less<>{},
664 [](const auto& a) { return a.first; });
665 if (it != sd->second.end() and it->first == id)
666 itg.first->second.emplace_back(id, k, it->second, active_coeffs);
667 }
668
669 if (integral->needs_facet_permutations)
670 needs_facet_permutations = true;
671 }
672 }
673
674 std::map<IntegralType,
675 std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>>
676 sd;
677 for (auto& [itg, data] : subdomains)
678 {
679 std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>> x;
680 for (auto& [id, idx] : data)
681 x.emplace_back(id, std::vector(idx.data(), idx.data() + idx.size()));
682 sd.insert({itg, std::move(x)});
683 }
684
685 return Form<T, U>(spaces, integrals, coefficients, constants,
686 needs_facet_permutations, entity_maps, mesh);
687}
688
702template <dolfinx::scalar T, std::floating_point U = scalar_value_type_t<T>>
704 const ufcx_form& ufcx_form,
705 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
706 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
707 coefficients,
708 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
709 const std::map<
711 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
712 subdomains,
713 const std::map<std::shared_ptr<const mesh::Mesh<U>>,
714 std::span<const std::int32_t>>& entity_maps,
715 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
716{
717 // Place coefficients in appropriate order
718 std::vector<std::shared_ptr<const Function<T, U>>> coeff_map;
719 for (const std::string& name : get_coefficient_names(ufcx_form))
720 {
721 if (auto it = coefficients.find(name); it != coefficients.end())
722 coeff_map.push_back(it->second);
723 else
724 {
725 throw std::runtime_error("Form coefficient \"" + name
726 + "\" not provided.");
727 }
728 }
729
730 // Place constants in appropriate order
731
732 std::vector<std::shared_ptr<const Constant<T>>> const_map;
733 for (const std::string& name : get_constant_names(ufcx_form))
734 {
735 if (auto it = constants.find(name); it != constants.end())
736 const_map.push_back(it->second);
737 else
738 throw std::runtime_error("Form constant \"" + name + "\" not provided.");
739 }
740
741 return create_form_factory(ufcx_form, spaces, coeff_map, const_map,
742 subdomains, entity_maps, mesh);
743}
744
762template <dolfinx::scalar T, std::floating_point U = scalar_value_type_t<T>>
764 ufcx_form* (*fptr)(),
765 const std::vector<std::shared_ptr<const FunctionSpace<U>>>& spaces,
766 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
767 coefficients,
768 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
769 const std::map<
771 std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>&
772 subdomains,
773 const std::map<std::shared_ptr<const mesh::Mesh<U>>,
774 std::span<const std::int32_t>>& entity_maps,
775 std::shared_ptr<const mesh::Mesh<U>> mesh = nullptr)
776{
777 ufcx_form* form = fptr();
778 Form<T, U> L = create_form<T, U>(*form, spaces, coefficients, constants,
779 subdomains, entity_maps, mesh);
780 std::free(form);
781 return L;
782}
783
795template <std::floating_point T>
797 std::shared_ptr<mesh::Mesh<T>> mesh, const basix::FiniteElement<T>& e,
798 const std::vector<std::size_t>& value_shape = {},
799 std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
800 reorder_fn
801 = nullptr)
802{
803 if (!e.value_shape().empty() and !value_shape.empty())
804 {
805 throw std::runtime_error(
806 "Cannot specify value shape for non-scalar base element.");
807 }
808
809 if (mesh::cell_type_from_basix_type(e.cell_type())
810 != mesh->topology()->cell_type())
811 throw std::runtime_error("Cell type of element and mesh must match.");
812
813 std::size_t bs = value_shape.empty()
814 ? 1
815 : std::accumulate(value_shape.begin(), value_shape.end(),
816 1, std::multiplies{});
817
818 // Create a DOLFINx element
819 auto _e = std::make_shared<const FiniteElement<T>>(e, bs);
820 assert(_e);
821
822 const std::vector<std::size_t> _value_shape
823 = (value_shape.empty() and !e.value_shape().empty())
824 ? fem::compute_value_shape(_e, mesh->topology()->dim(),
825 mesh->geometry().dim())
826 : value_shape;
827
828 // Create UFC subdofmaps and compute offset
829 const int num_sub_elements = _e->num_sub_elements();
830 std::vector<ElementDofLayout> sub_doflayout;
831 sub_doflayout.reserve(num_sub_elements);
832 for (int i = 0; i < num_sub_elements; ++i)
833 {
834 auto sub_element = _e->extract_sub_element({i});
835 std::vector<int> parent_map_sub(sub_element->space_dimension());
836 for (std::size_t j = 0; j < parent_map_sub.size(); ++j)
837 parent_map_sub[j] = i + _e->block_size() * j;
838 sub_doflayout.emplace_back(1, e.entity_dofs(), e.entity_closure_dofs(),
839 parent_map_sub, std::vector<ElementDofLayout>());
840 }
841
842 // Create a dofmap
843 ElementDofLayout layout(_e->block_size(), e.entity_dofs(),
844 e.entity_closure_dofs(), {}, sub_doflayout);
845 std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv
846 = nullptr;
847 if (_e->needs_dof_permutations())
848 permute_inv = _e->dof_permutation_fn(true, true);
849 assert(mesh);
850 assert(mesh->topology());
851 auto dofmap = std::make_shared<const DofMap>(create_dofmap(
852 mesh->comm(), layout, *mesh->topology(), permute_inv, reorder_fn));
853 return FunctionSpace(mesh, _e, dofmap, _value_shape);
854}
855
857namespace impl
858{
860template <dolfinx::scalar T, std::floating_point U>
861std::span<const std::uint32_t>
862get_cell_orientation_info(const Function<T, U>& coefficient)
863{
864 std::span<const std::uint32_t> cell_info;
865 auto element = coefficient.function_space()->element();
866 if (element->needs_dof_transformations())
867 {
868 auto mesh = coefficient.function_space()->mesh();
869 mesh->topology_mutable()->create_entity_permutations();
870 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
871 }
872
873 return cell_info;
874}
875
877template <int _bs, dolfinx::scalar T>
878void pack(std::span<T> coeffs, std::int32_t cell, int bs, std::span<const T> v,
879 std::span<const std::uint32_t> cell_info, const DofMap& dofmap,
880 auto transform)
881{
882 auto dofs = dofmap.cell_dofs(cell);
883 for (std::size_t i = 0; i < dofs.size(); ++i)
884 {
885 if constexpr (_bs < 0)
886 {
887 const int pos_c = bs * i;
888 const int pos_v = bs * dofs[i];
889 for (int k = 0; k < bs; ++k)
890 coeffs[pos_c + k] = v[pos_v + k];
891 }
892 else
893 {
894 assert(_bs == bs);
895 const int pos_c = _bs * i;
896 const int pos_v = _bs * dofs[i];
897 for (int k = 0; k < _bs; ++k)
898 coeffs[pos_c + k] = v[pos_v + k];
899 }
900 }
901
902 transform(coeffs, cell_info, cell, 1);
903}
904
907template <typename F>
908concept FetchCells = requires(F&& f, std::span<const std::int32_t> v) {
909 requires std::invocable<F, std::span<const std::int32_t>>;
910 { f(v) } -> std::convertible_to<std::int32_t>;
911};
912
926template <dolfinx::scalar T, std::floating_point U>
927void pack_coefficient_entity(std::span<T> c, int cstride,
928 const Function<T, U>& u,
929 std::span<const std::uint32_t> cell_info,
930 std::span<const std::int32_t> entities,
931 std::size_t estride, FetchCells auto&& fetch_cells,
932 std::int32_t offset)
933{
934 // Read data from coefficient Function u
935 std::span<const T> v = u.x()->array();
936 const DofMap& dofmap = *u.function_space()->dofmap();
937 auto element = u.function_space()->element();
938 assert(element);
939 int space_dim = element->space_dimension();
940
941 // Transformation from conforming degrees-of-freedom to reference
942 // degrees-of-freedom
943 auto transformation
944 = element->template dof_transformation_fn<T>(doftransform::transpose);
945 const int bs = dofmap.bs();
946 switch (bs)
947 {
948 case 1:
949 for (std::size_t e = 0; e < entities.size(); e += estride)
950 {
951 auto entity = entities.subspan(e, estride);
952 if (std::int32_t cell = fetch_cells(entity); cell >= 0)
953 {
954 auto cell_coeff
955 = c.subspan((e / estride) * cstride + offset, space_dim);
956 pack<1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
957 }
958 }
959 break;
960 case 2:
961 for (std::size_t e = 0; e < entities.size(); e += estride)
962 {
963 auto entity = entities.subspan(e, estride);
964 if (std::int32_t cell = fetch_cells(entity); cell >= 0)
965 {
966 auto cell_coeff
967 = c.subspan((e / estride) * cstride + offset, space_dim);
968 pack<2>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
969 }
970 }
971 break;
972 case 3:
973 for (std::size_t e = 0; e < entities.size(); e += estride)
974 {
975 auto entity = entities.subspan(e, estride);
976 if (std::int32_t cell = fetch_cells(entity); cell >= 0)
977 {
978 auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim);
979 pack<3>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
980 }
981 }
982 break;
983 default:
984 for (std::size_t e = 0; e < entities.size(); e += estride)
985 {
986 auto entity = entities.subspan(e, estride);
987 if (std::int32_t cell = fetch_cells(entity); cell >= 0)
988 {
989 auto cell_coeff
990 = c.subspan((e / estride) * cstride + offset, space_dim);
991 pack<-1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
992 }
993 }
994 break;
995 }
996}
997
998} // namespace impl
999
1006template <dolfinx::scalar T, std::floating_point U>
1007std::pair<std::vector<T>, int>
1009 int id)
1010{
1011 // Get form coefficient offsets and dofmaps
1012 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients
1013 = form.coefficients();
1014 const std::vector<int> offsets = form.coefficient_offsets();
1015
1016 std::size_t num_entities = 0;
1017 int cstride = 0;
1018 if (!coefficients.empty())
1019 {
1020 cstride = offsets.back();
1021 num_entities = form.domain(integral_type, id).size();
1022 if (integral_type == IntegralType::exterior_facet
1023 or integral_type == IntegralType::interior_facet)
1024 {
1025 num_entities /= 2;
1026 }
1027 }
1028
1029 return {std::vector<T>(num_entities * cstride), cstride};
1030}
1031
1036template <dolfinx::scalar T, std::floating_point U>
1037std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>
1039{
1040 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>> coeffs;
1041 for (auto integral_type : form.integral_types())
1042 {
1043 for (int id : form.integral_ids(integral_type))
1044 {
1045 coeffs.emplace_hint(
1046 coeffs.end(), std::pair(integral_type, id),
1047 allocate_coefficient_storage(form, integral_type, id));
1048 }
1049 }
1050
1051 return coeffs;
1052}
1053
1061template <dolfinx::scalar T, std::floating_point U>
1062void pack_coefficients(const Form<T, U>& form, IntegralType integral_type,
1063 int id, std::span<T> c, int cstride)
1064{
1065 // Get form coefficient offsets and dofmaps
1066 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients
1067 = form.coefficients();
1068 const std::vector<int> offsets = form.coefficient_offsets();
1069
1070 // Indicator for packing coefficients
1071 std::vector<int> active_coefficient(coefficients.size(), 0);
1072 if (!coefficients.empty())
1073 {
1074 switch (integral_type)
1075 {
1076 case IntegralType::cell:
1077 {
1078 // Get indicator for all coefficients that are active in cell
1079 // integrals
1080 for (std::size_t i = 0; i < form.num_integrals(IntegralType::cell); ++i)
1081 {
1082 for (auto idx : form.active_coeffs(IntegralType::cell, i))
1083 active_coefficient[idx] = 1;
1084 }
1085
1086 // Iterate over coefficients
1087 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
1088 {
1089 if (!active_coefficient[coeff])
1090 continue;
1091
1092 // Get coefficient mesh
1093 auto mesh = coefficients[coeff]->function_space()->mesh();
1094 assert(mesh);
1095
1096 // Other integrals in the form might have coefficients defined over
1097 // entities of codim > 0, which don't make sense for cell integrals, so
1098 // don't pack them.
1099 if (const int codim
1100 = form.mesh()->topology()->dim() - mesh->topology()->dim();
1101 codim > 0)
1102 {
1103 throw std::runtime_error("Should not be packing coefficients with "
1104 "codim>0 in a cell integral");
1105 }
1106
1107 std::vector<std::int32_t> cells
1108 = form.domain(IntegralType::cell, id, *mesh);
1109 std::span<const std::uint32_t> cell_info
1110 = impl::get_cell_orientation_info(*coefficients[coeff]);
1111 impl::pack_coefficient_entity(
1112 c, cstride, *coefficients[coeff], cell_info, cells, 1,
1113 [](auto entity) { return entity.front(); }, offsets[coeff]);
1114 }
1115 break;
1116 }
1118 {
1119 // Get indicator for all coefficients that are active in exterior
1120 // facet integrals
1121 for (std::size_t i = 0;
1123 {
1124 for (auto idx : form.active_coeffs(IntegralType::exterior_facet, i))
1125 active_coefficient[idx] = 1;
1126 }
1127
1128 // Iterate over coefficients
1129 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
1130 {
1131 if (!active_coefficient[coeff])
1132 continue;
1133
1134 auto mesh = coefficients[coeff]->function_space()->mesh();
1135 std::vector<std::int32_t> facets
1136 = form.domain(IntegralType::exterior_facet, id, *mesh);
1137 std::span<const std::uint32_t> cell_info
1138 = impl::get_cell_orientation_info(*coefficients[coeff]);
1139 impl::pack_coefficient_entity(
1140 c, cstride, *coefficients[coeff], cell_info, facets, 2,
1141 [](auto entity) { return entity.front(); }, offsets[coeff]);
1142 }
1143 break;
1144 }
1146 {
1147 // Get indicator for all coefficients that are active in interior
1148 // facet integrals
1149 for (std::size_t i = 0;
1151 {
1152 for (auto idx : form.active_coeffs(IntegralType::interior_facet, i))
1153 active_coefficient[idx] = 1;
1154 }
1155
1156 // Iterate over coefficients
1157 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
1158 {
1159 if (!active_coefficient[coeff])
1160 continue;
1161
1162 auto mesh = coefficients[coeff]->function_space()->mesh();
1163 std::vector<std::int32_t> facets
1164 = form.domain(IntegralType::interior_facet, id, *mesh);
1165
1166 std::span<const std::uint32_t> cell_info
1167 = impl::get_cell_orientation_info(*coefficients[coeff]);
1168
1169 // Pack coefficient ['+']
1170 impl::pack_coefficient_entity(
1171 c, 2 * cstride, *coefficients[coeff], cell_info, facets, 4,
1172 [](auto entity) { return entity[0]; }, 2 * offsets[coeff]);
1173
1174 // Pack coefficient ['-']
1175 impl::pack_coefficient_entity(
1176 c, 2 * cstride, *coefficients[coeff], cell_info, facets, 4,
1177 [](auto entity) { return entity[2]; },
1178 offsets[coeff] + offsets[coeff + 1]);
1179 }
1180 break;
1181 }
1182 default:
1183 throw std::runtime_error(
1184 "Could not pack coefficient. Integral type not supported.");
1185 }
1186 }
1187}
1188
1190template <dolfinx::scalar T, std::floating_point U = scalar_value_type_t<T>>
1192 const ufcx_expression& e,
1193 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients,
1194 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
1195 std::shared_ptr<const FunctionSpace<U>> argument_function_space = nullptr)
1196{
1197 if (e.rank > 0 and !argument_function_space)
1198 {
1199 throw std::runtime_error("Expression has Argument but no Argument "
1200 "function space was provided.");
1201 }
1202
1203 std::vector<U> X(e.points, e.points + e.num_points * e.entity_dimension);
1204 std::array<std::size_t, 2> Xshape
1205 = {static_cast<std::size_t>(e.num_points),
1206 static_cast<std::size_t>(e.entity_dimension)};
1207 std::vector<int> value_shape(e.value_shape, e.value_shape + e.num_components);
1208 std::function<void(T*, const T*, const T*,
1209 const typename scalar_value_type<T>::value_type*,
1210 const int*, const std::uint8_t*)>
1211 tabulate_tensor = nullptr;
1212 if constexpr (std::is_same_v<T, float>)
1213 tabulate_tensor = e.tabulate_tensor_float32;
1214#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
1215 else if constexpr (std::is_same_v<T, std::complex<float>>)
1216 {
1217 tabulate_tensor = reinterpret_cast<void (*)(
1218 T*, const T*, const T*,
1219 const typename scalar_value_type<T>::value_type*, const int*,
1220 const unsigned char*)>(e.tabulate_tensor_complex64);
1221 }
1222#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
1223 else if constexpr (std::is_same_v<T, double>)
1224 tabulate_tensor = e.tabulate_tensor_float64;
1225#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
1226 else if constexpr (std::is_same_v<T, std::complex<double>>)
1227 {
1228 tabulate_tensor = reinterpret_cast<void (*)(
1229 T*, const T*, const T*,
1230 const typename scalar_value_type<T>::value_type*, const int*,
1231 const unsigned char*)>(e.tabulate_tensor_complex128);
1232 }
1233#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS
1234 else
1235 throw std::runtime_error("Type not supported.");
1236
1237 assert(tabulate_tensor);
1238 return Expression(coefficients, constants, std::span<const U>(X), Xshape,
1239 tabulate_tensor, value_shape, argument_function_space);
1240}
1241
1244template <dolfinx::scalar T, std::floating_point U = scalar_value_type_t<T>>
1246 const ufcx_expression& e,
1247 const std::map<std::string, std::shared_ptr<const Function<T, U>>>&
1248 coefficients,
1249 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
1250 std::shared_ptr<const FunctionSpace<U>> argument_function_space = nullptr)
1251{
1252 // Place coefficients in appropriate order
1253 std::vector<std::shared_ptr<const Function<T, U>>> coeff_map;
1254 std::vector<std::string> coefficient_names;
1255 for (int i = 0; i < e.num_coefficients; ++i)
1256 coefficient_names.push_back(e.coefficient_names[i]);
1257
1258 for (const std::string& name : coefficient_names)
1259 {
1260 if (auto it = coefficients.find(name); it != coefficients.end())
1261 coeff_map.push_back(it->second);
1262 else
1263 {
1264 throw std::runtime_error("Expression coefficient \"" + name
1265 + "\" not provided.");
1266 }
1267 }
1268
1269 // Place constants in appropriate order
1270 std::vector<std::shared_ptr<const Constant<T>>> const_map;
1271 std::vector<std::string> constant_names;
1272 for (int i = 0; i < e.num_constants; ++i)
1273 constant_names.push_back(e.constant_names[i]);
1274
1275 for (const std::string& name : constant_names)
1276 {
1277 if (auto it = constants.find(name); it != constants.end())
1278 const_map.push_back(it->second);
1279 else
1280 {
1281 throw std::runtime_error("Expression constant \"" + name
1282 + "\" not provided.");
1283 }
1284 }
1285
1286 return create_expression(e, coeff_map, const_map, argument_function_space);
1287}
1288
1299template <dolfinx::scalar T, std::floating_point U>
1301 std::map<std::pair<IntegralType, int>,
1302 std::pair<std::vector<T>, int>>& coeffs)
1303{
1304 for (auto& [key, val] : coeffs)
1305 pack_coefficients<T>(form, key.first, key.second, val.first, val.second);
1306}
1307
1316template <dolfinx::scalar T, std::floating_point U>
1317std::pair<std::vector<T>, int>
1319 std::span<const std::int32_t> entities, std::size_t estride)
1320{
1321 // Get form coefficient offsets and dofmaps
1322 const std::vector<std::shared_ptr<const Function<T, U>>>& coeffs
1323 = e.coefficients();
1324 const std::vector<int> offsets = e.coefficient_offsets();
1325
1326 // Copy data into coefficient array
1327 const int cstride = offsets.back();
1328 std::vector<T> c(entities.size() / estride * offsets.back());
1329 if (!coeffs.empty())
1330 {
1331 // Iterate over coefficients
1332 for (std::size_t coeff = 0; coeff < coeffs.size(); ++coeff)
1333 {
1334 std::span<const std::uint32_t> cell_info
1335 = impl::get_cell_orientation_info(*coeffs[coeff]);
1336
1337 impl::pack_coefficient_entity(
1338 std::span(c), cstride, *coeffs[coeff], cell_info, entities, estride,
1339 [](auto entity) { return entity[0]; }, offsets[coeff]);
1340 }
1341 }
1342 return {std::move(c), cstride};
1343}
1344
1347template <typename U>
1348std::vector<typename U::scalar_type> pack_constants(const U& u)
1349{
1350 using T = typename U::scalar_type;
1351 const std::vector<std::shared_ptr<const Constant<T>>>& constants
1352 = u.constants();
1353
1354 // Calculate size of array needed to store packed constants
1355 std::int32_t size = std::accumulate(constants.cbegin(), constants.cend(), 0,
1356 [](std::int32_t sum, auto& constant)
1357 { return sum + constant->value.size(); });
1358
1359 // Pack constants
1360 std::vector<T> constant_values(size);
1361 std::int32_t offset = 0;
1362 for (auto& constant : constants)
1363 {
1364 const std::vector<T>& value = constant->value;
1365 std::ranges::copy(value, std::next(constant_values.begin(), offset));
1366 offset += value.size();
1367 }
1368
1369 return constant_values;
1370}
1371
1372} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
Definition utils.h:42
Definition Timer.h:32
double stop()
Definition Timer.cpp:38
Constant value which can be attached to a Form.
Definition Form.h:29
Degree-of-freedom map.
Definition DofMap.h:76
std::span< const std::int32_t > cell_dofs(std::int32_t c) const
Local-to-global mapping of dofs on a cell.
Definition DofMap.h:130
int bs() const noexcept
Return the block size for the dofmap.
Definition DofMap.cpp:173
Definition ElementDofLayout.h:30
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell.
Definition Function.h:32
Model of a finite element.
Definition FiniteElement.h:38
A representation of finite element variational forms.
Definition Form.h:139
const std::vector< std::shared_ptr< const Function< scalar_type, geometry_type > > > & coefficients() const
Access coefficients.
Definition Form.h:474
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell.
Definition Form.h:488
std::span< const std::int32_t > domain(IntegralType type, int i) const
Get the list of mesh entity indices for the ith integral (kernel) of a given type.
Definition Form.h:350
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
Extract common mesh for the form.
Definition Form.h:248
std::vector< int > integral_ids(IntegralType type) const
Get the IDs for integrals (kernels) for given integral type.
Definition Form.h:324
std::vector< int > active_coeffs(IntegralType type, std::size_t i) const
Indices of coefficients that are active for a given integral (kernel).
Definition Form.h:311
std::set< IntegralType > integral_types() const
Get types of integrals in the form.
Definition Form.h:281
int num_integrals(IntegralType type) const
Number of integrals on given domain type.
Definition Form.h:296
This class represents a finite element function space defined by a mesh, a finite element,...
Definition vtk_utils.h:32
Definition XDMFFile.h:29
Definition topologycomputation.h:24
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:44
Concepts for function that returns cell index.
Definition utils.h:908
void pack(std::span< T > coeffs, std::int32_t cell, int bs, std::span< const T > v, std::span< const std::uint32_t > cell_info, const DofMap &dofmap, auto transform)
Pack a single coefficient for a single cell.
Definition utils.h:878
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:63
void pack_coefficient_entity(std::span< T > c, int cstride, const Function< T, U > &u, std::span< const std::uint32_t > cell_info, std::span< const std::int32_t > entities, std::size_t estride, FetchCells auto &&fetch_cells, std::int32_t offset)
Pack a single coefficient for a set of active entities.
Definition utils.h:927
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_matrix_impl.h:26
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
std::pair< std::vector< T >, int > allocate_coefficient_storage(const Form< T, U > &form, IntegralType integral_type, int id)
Allocate storage for coefficients of a pair (integral_type, / id) from a Form.
Definition utils.h:1008
Form< T, U > create_form_factory(const ufcx_form &ufcx_form, 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:333
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
FunctionSpace(U mesh, V element, W dofmap, X value_shape) -> FunctionSpace< typename std::remove_cvref< typename U::element_type >::type::geometry_type::value_type >
Type deduction.
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:122
std::vector< std::int32_t > compute_integration_domains(IntegralType integral_type, const mesh::Topology &topology, std::span< const std::int32_t > entities, int dim)
Given an integral type and a set of entities, compute the entities that should be integrated over.
Definition utils.cpp:131
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u into a single array ready for assembly.
Definition utils.h:1348
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_function_space=nullptr)
Create Expression from UFC.
Definition utils.h:1191
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:241
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:703
IntegralType
Type of integral.
Definition Form.h:35
@ interior_facet
Interior facet.
@ exterior_facet
Exterior facet.
FunctionSpace< T > create_functionspace(std::shared_ptr< mesh::Mesh< T > > mesh, const basix::FiniteElement< T > &e, const std::vector< std::size_t > &value_shape={}, std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> reorder_fn=nullptr)
Create a function space from a Basix element.
Definition utils.h:796
la::SparsityPattern create_sparsity_pattern(const Form< T, U > &a)
Create a sparsity pattern for a given form.
Definition utils.h:147
std::vector< std::size_t > compute_value_shape(std::shared_ptr< const dolfinx::fem::FiniteElement< T > > element, std::size_t tdim, std::size_t gdim)
Compute the physical value shape of an element for a mesh.
Definition FunctionSpace.h:439
void pack_coefficients(const Form< T, U > &form, IntegralType integral_type, int id, std::span< T > c, int cstride)
Pack coefficients of a Form for a given integral type and domain id.
Definition utils.h:1062
CellType cell_type_from_basix_type(basix::cell::type celltype)
Get a cell type from a Basix cell type.
Definition cell_types.cpp:237
std::vector< std::int32_t > exterior_facet_indices(const Topology &topology)
Compute the indices of all exterior facets that are owned by the caller.
Definition utils.cpp:58