DOLFINx
DOLFINx C++ interface
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 "sparsitybuild.h"
17#include <array>
18#include <dolfinx/la/SparsityPattern.h>
19#include <dolfinx/mesh/cell_types.h>
20#include <functional>
21#include <memory>
22#include <set>
23#include <span>
24#include <stdexcept>
25#include <string>
26#include <type_traits>
27#include <ufcx.h>
28#include <utility>
29#include <vector>
30
33
34namespace basix
35{
36class FiniteElement;
37}
38
39namespace dolfinx::common
40{
41class IndexMap;
42}
43
44namespace dolfinx::mesh
45{
46class Mesh;
47class Topology;
48} // namespace dolfinx::mesh
49
50namespace dolfinx::fem
51{
52class FunctionSpace;
53
54namespace impl
55{
58template <typename T, typename = void>
59struct scalar_value_type
60{
62 typedef T value_type;
63};
65template <typename T>
66struct scalar_value_type<T, std::void_t<typename T::value_type>>
67{
68 typedef typename T::value_type value_type;
69};
71template <typename T>
72using scalar_value_type_t = typename scalar_value_type<T>::value_type;
73
74} // namespace impl
75
83template <typename T>
84std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace>, 2>>>
85extract_function_spaces(const std::vector<std::vector<const Form<T>*>>& a)
86{
87 std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace>, 2>>>
88 spaces(a.size(),
89 std::vector<std::array<std::shared_ptr<const FunctionSpace>, 2>>(
90 a[0].size()));
91 for (std::size_t i = 0; i < a.size(); ++i)
92 {
93 for (std::size_t j = 0; j < a[i].size(); ++j)
94 {
95 if (const Form<T>* form = a[i][j]; form)
96 spaces[i][j] = {form->function_spaces()[0], form->function_spaces()[1]};
97 }
98 }
99 return spaces;
100}
101
106 const mesh::Topology& topology,
107 const std::array<std::reference_wrapper<const DofMap>, 2>& dofmaps,
108 const std::set<IntegralType>& integrals);
109
115template <typename T>
117{
118 if (a.rank() != 2)
119 {
120 throw std::runtime_error(
121 "Cannot create sparsity pattern. Form is not a bilinear form");
122 }
123
124 // Get dof maps and mesh
125 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
126 *a.function_spaces().at(0)->dofmap(),
127 *a.function_spaces().at(1)->dofmap()};
128 std::shared_ptr mesh = a.mesh();
129 assert(mesh);
130
131 const std::set<IntegralType> types = a.integral_types();
132 if (types.find(IntegralType::interior_facet) != types.end()
133 or types.find(IntegralType::exterior_facet) != types.end())
134 {
135 // FIXME: cleanup these calls? Some of the happen internally again.
136 const int tdim = mesh->topology().dim();
137 mesh->topology_mutable().create_entities(tdim - 1);
138 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
139 }
140
141 common::Timer t0("Build sparsity");
142
143 // Get common::IndexMaps for each dimension
144 const std::array index_maps{dofmaps[0].get().index_map,
145 dofmaps[1].get().index_map};
146 const std::array bs
147 = {dofmaps[0].get().index_map_bs(), dofmaps[1].get().index_map_bs()};
148
149 // Create and build sparsity pattern
150 la::SparsityPattern pattern(mesh->comm(), index_maps, bs);
151 for (auto type : types)
152 {
153 std::vector<int> ids = a.integral_ids(type);
154 switch (type)
155 {
157 for (int id : ids)
158 {
159 const std::vector<std::int32_t>& cells = a.cell_domains(id);
160 sparsitybuild::cells(pattern, cells, {{dofmaps[0], dofmaps[1]}});
161 }
162 break;
164 for (int id : ids)
165 {
166 const std::vector<std::int32_t>& facets = a.interior_facet_domains(id);
167 sparsitybuild::interior_facets(pattern, facets,
168 {{dofmaps[0], dofmaps[1]}});
169 }
170 break;
172 for (int id : ids)
173 {
174 const std::vector<std::int32_t>& facets = a.exterior_facet_domains(id);
175 sparsitybuild::exterior_facets(pattern, facets,
176 {{dofmaps[0], dofmaps[1]}});
177 }
178 break;
179 default:
180 throw std::runtime_error("Unsupported integral type");
181 }
182 }
183
184 t0.stop();
185
186 return pattern;
187}
188
190ElementDofLayout create_element_dof_layout(const ufcx_dofmap& dofmap,
191 const mesh::CellType cell_type,
192 const std::vector<int>& parent_map
193 = {});
194
203DofMap
204create_dofmap(MPI_Comm comm, const ElementDofLayout& layout,
205 mesh::Topology& topology,
206 const std::function<std::vector<int>(
207 const graph::AdjacencyList<std::int32_t>&)>& reorder_fn,
208 const FiniteElement& element);
209
213std::vector<std::string> get_coefficient_names(const ufcx_form& ufcx_form);
214
218std::vector<std::string> get_constant_names(const ufcx_form& ufcx_form);
219
227template <typename T>
229 const ufcx_form& ufcx_form,
230 const std::vector<std::shared_ptr<const FunctionSpace>>& spaces,
231 const std::vector<std::shared_ptr<const Function<T>>>& coefficients,
232 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
233 const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
234 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
235{
236 if (ufcx_form.rank != (int)spaces.size())
237 throw std::runtime_error("Wrong number of argument spaces for Form.");
238 if (ufcx_form.num_coefficients != (int)coefficients.size())
239 {
240 throw std::runtime_error(
241 "Mismatch between number of expected and provided Form coefficients.");
242 }
243 if (ufcx_form.num_constants != (int)constants.size())
244 {
245 throw std::runtime_error(
246 "Mismatch between number of expected and provided Form constants.");
247 }
248
249 // Check argument function spaces
250#ifndef NDEBUG
251 for (std::size_t i = 0; i < spaces.size(); ++i)
252 {
253 assert(spaces[i]->element());
254 ufcx_finite_element* ufcx_element = ufcx_form.finite_elements[i];
255 assert(ufcx_element);
256 if (std::string(ufcx_element->signature)
257 != spaces[i]->element()->signature())
258 {
259 throw std::runtime_error(
260 "Cannot create form. Wrong type of function space for argument.");
261 }
262 }
263#endif
264
265 // Get list of integral IDs, and load tabulate tensor into memory for
266 // each
267 using kern = std::function<void(
268 T*, const T*, const T*,
269 const typename impl::scalar_value_type<T>::value_type*, const int*,
270 const std::uint8_t*)>;
271 std::map<IntegralType, std::pair<std::vector<std::pair<int, kern>>,
272 const mesh::MeshTags<int>*>>
273 integral_data;
274
275 bool needs_facet_permutations = false;
276
277 // Attach cell kernels
278 std::vector<int> cell_integral_ids(ufcx_form.integral_ids(cell),
279 ufcx_form.integral_ids(cell)
280 + ufcx_form.num_integrals(cell));
281 for (int i = 0; i < ufcx_form.num_integrals(cell); ++i)
282 {
283 ufcx_integral* integral = ufcx_form.integrals(cell)[i];
284 assert(integral);
285
286 kern k = nullptr;
287 if constexpr (std::is_same_v<T, float>)
288 k = integral->tabulate_tensor_float32;
289 else if constexpr (std::is_same_v<T, std::complex<float>>)
290 {
291 k = reinterpret_cast<void (*)(
292 T*, const T*, const T*,
293 const typename impl::scalar_value_type<T>::value_type*, const int*,
294 const unsigned char*)>(integral->tabulate_tensor_complex64);
295 }
296 else if constexpr (std::is_same_v<T, double>)
297 k = integral->tabulate_tensor_float64;
298 else if constexpr (std::is_same_v<T, std::complex<double>>)
299 {
300 k = reinterpret_cast<void (*)(
301 T*, const T*, const T*,
302 const typename impl::scalar_value_type<T>::value_type*, const int*,
303 const unsigned char*)>(integral->tabulate_tensor_complex128);
304 }
305 assert(k);
306
307 integral_data[IntegralType::cell].first.emplace_back(cell_integral_ids[i],
308 k);
309 if (integral->needs_facet_permutations)
310 needs_facet_permutations = true;
311 }
312
313 // Attach cell subdomain data
314 if (auto it = subdomains.find(IntegralType::cell);
315 it != subdomains.end() and !cell_integral_ids.empty())
316 {
317 integral_data[IntegralType::cell].second = it->second;
318 }
319
320 // FIXME: Can facets be handled better?
321
322 // Create facets, if required
323 if (ufcx_form.num_integrals(exterior_facet) > 0
324 or ufcx_form.num_integrals(interior_facet) > 0)
325 {
326 if (!spaces.empty())
327 {
328 auto mesh = spaces[0]->mesh();
329 const int tdim = mesh->topology().dim();
330 spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
331 }
332 }
333
334 // Attach exterior facet kernels
335 std::vector<int> exterior_facet_integral_ids(
336 ufcx_form.integral_ids(exterior_facet),
337 ufcx_form.integral_ids(exterior_facet)
338 + ufcx_form.num_integrals(exterior_facet));
339 for (int i = 0; i < ufcx_form.num_integrals(exterior_facet); ++i)
340 {
341 ufcx_integral* integral = ufcx_form.integrals(exterior_facet)[i];
342 assert(integral);
343
344 kern k = nullptr;
345 if constexpr (std::is_same_v<T, float>)
346 k = integral->tabulate_tensor_float32;
347 else if constexpr (std::is_same_v<T, std::complex<float>>)
348 {
349 k = reinterpret_cast<void (*)(
350 T*, const T*, const T*,
351 const typename impl::scalar_value_type<T>::value_type*, const int*,
352 const unsigned char*)>(integral->tabulate_tensor_complex64);
353 }
354 else if constexpr (std::is_same_v<T, double>)
355 k = integral->tabulate_tensor_float64;
356 else if constexpr (std::is_same_v<T, std::complex<double>>)
357 {
358 k = reinterpret_cast<void (*)(
359 T*, const T*, const T*,
360 const typename impl::scalar_value_type<T>::value_type*, const int*,
361 const unsigned char*)>(integral->tabulate_tensor_complex128);
362 }
363 assert(k);
364
365 integral_data[IntegralType::exterior_facet].first.emplace_back(
366 exterior_facet_integral_ids[i], k);
367 if (integral->needs_facet_permutations)
368 needs_facet_permutations = true;
369 }
370
371 // Attach exterior facet subdomain data
372 if (auto it = subdomains.find(IntegralType::exterior_facet);
373 it != subdomains.end() and !exterior_facet_integral_ids.empty())
374 {
375 integral_data[IntegralType::exterior_facet].second = it->second;
376 }
377
378 // Attach interior facet kernels
379 std::vector<int> interior_facet_integral_ids(
380 ufcx_form.integral_ids(interior_facet),
381 ufcx_form.integral_ids(interior_facet)
382 + ufcx_form.num_integrals(interior_facet));
383 for (int i = 0; i < ufcx_form.num_integrals(interior_facet); ++i)
384 {
385 ufcx_integral* integral = ufcx_form.integrals(interior_facet)[i];
386 assert(integral);
387
388 kern k = nullptr;
389 if constexpr (std::is_same_v<T, float>)
390 k = integral->tabulate_tensor_float32;
391 else if constexpr (std::is_same_v<T, std::complex<float>>)
392 {
393 k = reinterpret_cast<void (*)(
394 T*, const T*, const T*,
395 const typename impl::scalar_value_type<T>::value_type*, const int*,
396 const unsigned char*)>(integral->tabulate_tensor_complex64);
397 }
398 else if constexpr (std::is_same_v<T, double>)
399 k = integral->tabulate_tensor_float64;
400 else if constexpr (std::is_same_v<T, std::complex<double>>)
401 {
402 k = reinterpret_cast<void (*)(
403 T*, const T*, const T*,
404 const typename impl::scalar_value_type<T>::value_type*, const int*,
405 const unsigned char*)>(integral->tabulate_tensor_complex128);
406 }
407 assert(k);
408
409 integral_data[IntegralType::interior_facet].first.emplace_back(
410 interior_facet_integral_ids[i], k);
411 if (integral->needs_facet_permutations)
412 needs_facet_permutations = true;
413 }
414
415 // Attach interior facet subdomain data
416 if (auto it = subdomains.find(IntegralType::interior_facet);
417 it != subdomains.end() and !interior_facet_integral_ids.empty())
418 {
419 integral_data[IntegralType::interior_facet].second = it->second;
420 }
421
422 return Form<T>(spaces, integral_data, coefficients, constants,
423 needs_facet_permutations, mesh);
424}
425
435template <typename T>
437 const ufcx_form& ufcx_form,
438 const std::vector<std::shared_ptr<const FunctionSpace>>& spaces,
439 const std::map<std::string, std::shared_ptr<const Function<T>>>&
440 coefficients,
441 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
442 const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
443 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
444{
445 // Place coefficients in appropriate order
446 std::vector<std::shared_ptr<const Function<T>>> coeff_map;
447 for (const std::string& name : get_coefficient_names(ufcx_form))
448 {
449 if (auto it = coefficients.find(name); it != coefficients.end())
450 coeff_map.push_back(it->second);
451 else
452 {
453 throw std::runtime_error("Form coefficient \"" + name
454 + "\" not provided.");
455 }
456 }
457
458 // Place constants in appropriate order
459 std::vector<std::shared_ptr<const Constant<T>>> const_map;
460 for (const std::string& name : get_constant_names(ufcx_form))
461 {
462 if (auto it = constants.find(name); it != constants.end())
463 const_map.push_back(it->second);
464 else
465 throw std::runtime_error("Form constant \"" + name + "\" not provided.");
466 }
467
468 return create_form(ufcx_form, spaces, coeff_map, const_map, subdomains, mesh);
469}
470
482template <typename T>
484 ufcx_form* (*fptr)(),
485 const std::vector<std::shared_ptr<const FunctionSpace>>& spaces,
486 const std::map<std::string, std::shared_ptr<const Function<T>>>&
487 coefficients,
488 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
489 const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
490 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
491{
492 ufcx_form* form = fptr();
493 Form<T> L = create_form<T>(*form, spaces, coefficients, constants, subdomains,
494 mesh);
495 std::free(form);
496 return L;
497}
498
507FunctionSpace
508create_functionspace(const std::shared_ptr<mesh::Mesh>& mesh,
509 const basix::FiniteElement& e, int bs,
510 const std::function<std::vector<int>(
511 const graph::AdjacencyList<std::int32_t>&)>& reorder_fn
512 = nullptr);
513
526FunctionSpace create_functionspace(
527 ufcx_function_space* (*fptr)(const char*), const std::string& function_name,
528 const std::shared_ptr<mesh::Mesh>& mesh,
529 const std::function<
530 std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>& reorder_fn
531 = nullptr);
532
534namespace impl
535{
537template <typename T>
538std::span<const std::uint32_t> get_cell_orientation_info(
539 const std::vector<std::shared_ptr<const Function<T>>>& coefficients)
540{
541 bool needs_dof_transformations = false;
542 for (auto coeff : coefficients)
543 {
544 std::shared_ptr<const FiniteElement> element
545 = coeff->function_space()->element();
546 if (element->needs_dof_transformations())
547 {
548 needs_dof_transformations = true;
549 break;
550 }
551 }
552
553 std::span<const std::uint32_t> cell_info;
554 if (needs_dof_transformations)
555 {
556 auto mesh = coefficients.front()->function_space()->mesh();
557 mesh->topology_mutable().create_entity_permutations();
558 cell_info = std::span(mesh->topology().get_cell_permutation_info());
559 }
560
561 return cell_info;
562}
563
564// Pack a single coefficient for a single cell
565template <typename T, int _bs, typename Functor>
566static inline void pack(const std::span<T>& coeffs, std::int32_t cell, int bs,
567 const std::span<const T>& v,
568 const std::span<const std::uint32_t>& cell_info,
569 const DofMap& dofmap, Functor transform)
570{
571 auto dofs = dofmap.cell_dofs(cell);
572 for (std::size_t i = 0; i < dofs.size(); ++i)
573 {
574 if constexpr (_bs < 0)
575 {
576 const int pos_c = bs * i;
577 const int pos_v = bs * dofs[i];
578 for (int k = 0; k < bs; ++k)
579 coeffs[pos_c + k] = v[pos_v + k];
580 }
581 else
582 {
583 const int pos_c = _bs * i;
584 const int pos_v = _bs * dofs[i];
585 for (int k = 0; k < _bs; ++k)
586 coeffs[pos_c + k] = v[pos_v + k];
587 }
588 }
589
590 transform(coeffs, cell_info, cell, 1);
591}
592
607template <typename T, typename Functor>
608void pack_coefficient_entity(const std::span<T>& c, int cstride,
609 const Function<T>& u,
610 const std::span<const std::uint32_t>& cell_info,
611 const std::span<const std::int32_t>& entities,
612 std::size_t estride, Functor fetch_cells,
613 std::int32_t offset)
614{
615 // Read data from coefficient "u"
616 const std::span<const T>& v = u.x()->array();
617 const DofMap& dofmap = *u.function_space()->dofmap();
618 std::shared_ptr<const FiniteElement> element = u.function_space()->element();
619 int space_dim = element->space_dimension();
620 const auto transformation
621 = element->get_dof_transformation_function<T>(false, true);
622
623 const int bs = dofmap.bs();
624 switch (bs)
625 {
626 case 1:
627 for (std::size_t e = 0; e < entities.size(); e += estride)
628 {
629 auto entity = entities.subspan(e, estride);
630 std::int32_t cell = fetch_cells(entity);
631 auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim);
632 pack<T, 1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
633 }
634 break;
635 case 2:
636 for (std::size_t e = 0; e < entities.size(); e += estride)
637 {
638 auto entity = entities.subspan(e, estride);
639 std::int32_t cell = fetch_cells(entity);
640 auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim);
641 pack<T, 2>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
642 }
643 break;
644 case 3:
645 for (std::size_t e = 0; e < entities.size(); e += estride)
646 {
647 auto entity = entities.subspan(e, estride);
648 std::int32_t cell = fetch_cells(entity);
649 auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim);
650 pack<T, 3>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
651 }
652 break;
653 default:
654 for (std::size_t e = 0; e < entities.size(); e += estride)
655 {
656 auto entity = entities.subspan(e, estride);
657 std::int32_t cell = fetch_cells(entity);
658 auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim);
659 pack<T, -1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
660 }
661 break;
662 }
663}
664
665} // namespace impl
666
673template <typename T>
674std::pair<std::vector<T>, int>
676 int id)
677{
678 // Get form coefficient offsets and dofmaps
679 const std::vector<std::shared_ptr<const Function<T>>>& coefficients
680 = form.coefficients();
681 const std::vector<int> offsets = form.coefficient_offsets();
682
683 std::size_t num_entities = 0;
684 int cstride = 0;
685 if (!coefficients.empty())
686 {
687 cstride = offsets.back();
688 switch (integral_type)
689 {
691 num_entities = form.cell_domains(id).size();
692 break;
694 num_entities = form.exterior_facet_domains(id).size() / 2;
695 break;
697 num_entities = form.interior_facet_domains(id).size() / 2;
698 break;
699 default:
700 throw std::runtime_error(
701 "Could not allocate coefficient data. Integral type not supported.");
702 }
703 }
704
705 return {std::vector<T>(num_entities * cstride), cstride};
706}
707
712template <typename T>
713std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>
715{
716 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>> coeffs;
717 for (auto integral_type : form.integral_types())
718 {
719 for (int id : form.integral_ids(integral_type))
720 {
721 coeffs.emplace_hint(
722 coeffs.end(), std::pair(integral_type, id),
723 allocate_coefficient_storage(form, integral_type, id));
724 }
725 }
726
727 return coeffs;
728}
729
737template <typename T>
738void pack_coefficients(const Form<T>& form, IntegralType integral_type, int id,
739 const std::span<T>& c, int cstride)
740{
741 // Get form coefficient offsets and dofmaps
742 const std::vector<std::shared_ptr<const Function<T>>>& coefficients
743 = form.coefficients();
744 const std::vector<int> offsets = form.coefficient_offsets();
745
746 if (!coefficients.empty())
747 {
748 std::span<const std::uint32_t> cell_info
749 = impl::get_cell_orientation_info(coefficients);
750
751 switch (integral_type)
752 {
754 {
755 auto fetch_cell = [](auto entity) { return entity.front(); };
756 const std::vector<std::int32_t>& cells = form.cell_domains(id);
757 // Iterate over coefficients
758 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
759 {
760 impl::pack_coefficient_entity(c, cstride, *coefficients[coeff],
761 cell_info, cells, 1, fetch_cell,
762 offsets[coeff]);
763 }
764 break;
765 }
767 {
768 const std::vector<std::int32_t>& facets = form.exterior_facet_domains(id);
769
770 // Create lambda function fetching cell index from exterior facet entity
771 auto fetch_cell = [](auto& entity) { return entity.front(); };
772
773 // Iterate over coefficients
774 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
775 {
776 impl::pack_coefficient_entity(c, cstride, *coefficients[coeff],
777 cell_info, facets, 2, fetch_cell,
778 offsets[coeff]);
779 }
780
781 break;
782 }
784 {
785 const std::vector<std::int32_t>& facets = form.interior_facet_domains(id);
786 // Lambda functions to fetch cell index from interior facet entity
787 auto fetch_cell0 = [](auto& entity) { return entity[0]; };
788 auto fetch_cell1 = [](auto& entity) { return entity[2]; };
789
790 // Iterate over coefficients
791 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
792 {
793 // Pack coefficient ['+']
794 impl::pack_coefficient_entity(c, 2 * cstride, *coefficients[coeff],
795 cell_info, facets, 4, fetch_cell0,
796 2 * offsets[coeff]);
797 // Pack coefficient ['-']
798 impl::pack_coefficient_entity(c, 2 * cstride, *coefficients[coeff],
799 cell_info, facets, 4, fetch_cell1,
800 offsets[coeff] + offsets[coeff + 1]);
801 }
802 break;
803 }
804 default:
805 throw std::runtime_error(
806 "Could not pack coefficient. Integral type not supported.");
807 }
808 }
809}
810
812template <typename T>
814 const ufcx_expression& expression,
815 const std::vector<std::shared_ptr<const Function<T>>>& coefficients,
816 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
817 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr,
818 const std::shared_ptr<const FunctionSpace>& argument_function_space
819 = nullptr)
820{
821 if (expression.rank > 0 and !argument_function_space)
822 {
823 throw std::runtime_error(
824 "Expression has Argument but no Argument function space was provided.");
825 }
826
827 const std::size_t size
828 = expression.num_points * expression.topological_dimension;
829 std::span<const double> X(expression.points, size);
830 std::array<std::size_t, 2> Xshape
831 = {static_cast<std::size_t>(expression.num_points),
832 static_cast<std::size_t>(expression.topological_dimension)};
833
834 std::vector<int> value_shape;
835 for (int i = 0; i < expression.num_components; ++i)
836 value_shape.push_back(expression.value_shape[i]);
837
838 std::function<void(T*, const T*, const T*,
839 const typename impl::scalar_value_type<T>::value_type*,
840 const int*, const std::uint8_t*)>
841 tabulate_tensor = nullptr;
842 if constexpr (std::is_same_v<T, float>)
843 tabulate_tensor = expression.tabulate_tensor_float32;
844 else if constexpr (std::is_same_v<T, std::complex<float>>)
845 {
846 tabulate_tensor = reinterpret_cast<void (*)(
847 T*, const T*, const T*,
848 const typename impl::scalar_value_type<T>::value_type*, const int*,
849 const unsigned char*)>(expression.tabulate_tensor_complex64);
850 }
851 else if constexpr (std::is_same_v<T, double>)
852 tabulate_tensor = expression.tabulate_tensor_float64;
853 else if constexpr (std::is_same_v<T, std::complex<double>>)
854 {
855 tabulate_tensor = reinterpret_cast<void (*)(
856 T*, const T*, const T*,
857 const typename impl::scalar_value_type<T>::value_type*, const int*,
858 const unsigned char*)>(expression.tabulate_tensor_complex128);
859 }
860 else
861 {
862 throw std::runtime_error("Type not supported.");
863 }
864 assert(tabulate_tensor);
865
866 return Expression(coefficients, constants, X, Xshape, tabulate_tensor,
867 value_shape, mesh, argument_function_space);
868}
869
872template <typename T>
874 const ufcx_expression& expression,
875 const std::map<std::string, std::shared_ptr<const Function<T>>>&
876 coefficients,
877 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
878 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr,
879 const std::shared_ptr<const FunctionSpace>& argument_function_space
880 = nullptr)
881{
882 // Place coefficients in appropriate order
883 std::vector<std::shared_ptr<const Function<T>>> coeff_map;
884 std::vector<std::string> coefficient_names;
885 for (int i = 0; i < expression.num_coefficients; ++i)
886 coefficient_names.push_back(expression.coefficient_names[i]);
887
888 for (const std::string& name : coefficient_names)
889 {
890 if (auto it = coefficients.find(name); it != coefficients.end())
891 coeff_map.push_back(it->second);
892 else
893 {
894 throw std::runtime_error("Expression coefficient \"" + name
895 + "\" not provided.");
896 }
897 }
898
899 // Place constants in appropriate order
900 std::vector<std::shared_ptr<const Constant<T>>> const_map;
901 std::vector<std::string> constant_names;
902 for (int i = 0; i < expression.num_constants; ++i)
903 constant_names.push_back(expression.constant_names[i]);
904
905 for (const std::string& name : constant_names)
906 {
907 if (auto it = constants.find(name); it != constants.end())
908 const_map.push_back(it->second);
909 else
910 {
911 throw std::runtime_error("Expression constant \"" + name
912 + "\" not provided.");
913 }
914 }
915
916 return create_expression(expression, coeff_map, const_map, mesh,
917 argument_function_space);
918}
919
925template <typename T>
926void pack_coefficients(const Form<T>& form,
927 std::map<std::pair<IntegralType, int>,
928 std::pair<std::vector<T>, int>>& coeffs)
929{
930 for (auto& [key, val] : coeffs)
931 pack_coefficients<T>(form, key.first, key.second, val.first, val.second);
932}
933
940template <typename T>
941std::pair<std::vector<T>, int>
943 const std::span<const std::int32_t>& cells)
944{
945 // Get form coefficient offsets and dofmaps
946 const std::vector<std::shared_ptr<const Function<T>>>& coefficients
947 = u.coefficients();
948 const std::vector<int> offsets = u.coefficient_offsets();
949
950 // Copy data into coefficient array
951 const int cstride = offsets.back();
952 std::vector<T> c(cells.size() * offsets.back());
953 if (!coefficients.empty())
954 {
955 std::span<const std::uint32_t> cell_info
956 = impl::get_cell_orientation_info(coefficients);
957
958 // Iterate over coefficients
959 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
960 impl::pack_coefficient_entity(
961 std::span(c), cstride, *coefficients[coeff], cell_info, cells, 1,
962 [](auto entity) { return entity[0]; }, offsets[coeff]);
963 }
964 return {std::move(c), cstride};
965}
966
969template <typename U>
970std::vector<typename U::scalar_type> pack_constants(const U& u)
971{
972 using T = typename U::scalar_type;
973 const std::vector<std::shared_ptr<const Constant<T>>>& constants
974 = u.constants();
975
976 // Calculate size of array needed to store packed constants
977 std::int32_t size = std::accumulate(constants.cbegin(), constants.cend(), 0,
978 [](std::int32_t sum, auto& constant)
979 { return sum + constant->value.size(); });
980
981 // Pack constants
982 std::vector<T> constant_values(size);
983 std::int32_t offset = 0;
984 for (auto& constant : constants)
985 {
986 const std::vector<T>& value = constant->value;
987 std::copy(value.cbegin(), value.cend(),
988 std::next(constant_values.begin(), offset));
989 offset += value.size();
990 }
991
992 return constant_values;
993}
994
995} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
A timer can be used for timing tasks. The basic usage is.
Definition: Timer.h:31
double stop()
Stop timer, return wall time elapsed and store timing data into logger.
Definition: Timer.cpp:45
Constant value which can be attached to a Form. Constants may be scalar (rank 0), vector (rank 1),...
Definition: Constant.h:20
Degree-of-freedom map.
Definition: DofMap.h:71
int bs() const noexcept
Return the block size for the dofmap.
Definition: DofMap.cpp:183
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell....
Definition: Expression.h:34
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell. Used to pack data for multiple coefficients in...
Definition: Expression.h:122
const std::vector< std::shared_ptr< const Function< T > > > & coefficients() const
Get coefficients.
Definition: Expression.h:105
A representation of finite element variational forms.
Definition: Form.h:63
int rank() const
Rank of the form (bilinear form = 2, linear form = 1, functional = 0, etc)
Definition: Form.h:172
const std::vector< std::int32_t > & cell_domains(int i) const
Get the list of cell indices for the ith integral (kernel) for the cell domain type.
Definition: Form.h:280
const std::vector< std::int32_t > & exterior_facet_domains(int i) const
Get the list of (cell_index, local_facet_index) pairs for the ith integral (kernel) for the exterior ...
Definition: Form.h:293
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell. Used to pack data for multiple coefficients in...
Definition: Form.h:330
const std::vector< std::shared_ptr< const FunctionSpace > > & function_spaces() const
Return function spaces for all arguments.
Definition: Form.h:181
const std::vector< std::shared_ptr< const Function< T > > > & coefficients() const
Access coefficients.
Definition: Form.h:317
std::vector< int > integral_ids(IntegralType type) const
Get the IDs for integrals (kernels) for given integral type. The IDs correspond to the domain IDs whi...
Definition: Form.h:248
std::set< IntegralType > integral_types() const
Get types of integrals in the form.
Definition: Form.h:211
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh for the form.
Definition: Form.h:176
const std::vector< std::int32_t > & interior_facet_domains(int i) const
Get the list of (cell_index_0, local_facet_index_0, cell_index_1, local_facet_index_1) quadruplets fo...
Definition: Form.h:308
This class represents a function in a finite element function space , given by.
Definition: Function.h:45
std::shared_ptr< const FunctionSpace > function_space() const
Access the function space.
Definition: Function.h:136
std::shared_ptr< const la::Vector< T > > x() const
Underlying vector.
Definition: Function.h:142
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:26
This class provides a sparsity pattern data structure that can be used to initialize sparse matrices....
Definition: SparsityPattern.h:34
MeshTags associate values with mesh entities.
Definition: MeshTags.h:36
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition: Topology.h:44
void pack_coefficient_entity(const std::span< T > &c, int cstride, const Function< T > &u, const std::span< const std::uint32_t > &cell_info, const std::span< const std::int32_t > &entities, std::size_t estride, Functor fetch_cells, std::int32_t offset)
Pack a single coefficient for a set of active entities.
Definition: utils.h:608
Miscellaneous classes, functions and types.
void interior_facets(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const std::reference_wrapper< const DofMap >, 2 > &dofmaps)
Iterate over interior facets and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:43
void cells(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const std::reference_wrapper< const DofMap >, 2 > &dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:18
void exterior_facets(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const std::reference_wrapper< const DofMap >, 2 > &dofmaps)
Iterate over exterior facets and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:112
Finite element method functionality.
Definition: assemble_matrix_impl.h:25
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:196
Expression< T > create_expression(const ufcx_expression &expression, const std::vector< std::shared_ptr< const Function< T > > > &coefficients, const std::vector< std::shared_ptr< const Constant< T > > > &constants, const std::shared_ptr< const mesh::Mesh > &mesh=nullptr, const std::shared_ptr< const FunctionSpace > &argument_function_space=nullptr)
Create Expression from UFC.
Definition: utils.h:813
std::vector< std::string > get_coefficient_names(const ufcx_form &ufcx_form)
Get the name of each coefficient in a UFC form.
Definition: utils.cpp:187
std::vector< std::vector< std::array< std::shared_ptr< const FunctionSpace >, 2 > > > extract_function_spaces(const std::vector< std::vector< const Form< T > * > > &a)
Extract test (0) and trial (1) function spaces pairs for each bilinear form for a rectangular array o...
Definition: utils.h:85
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:970
Form< T > create_form(const ufcx_form &ufcx_form, const std::vector< std::shared_ptr< const FunctionSpace > > &spaces, const std::vector< std::shared_ptr< const Function< T > > > &coefficients, const std::vector< std::shared_ptr< const Constant< T > > > &constants, const std::map< IntegralType, const mesh::MeshTags< int > * > &subdomains, const std::shared_ptr< const mesh::Mesh > &mesh=nullptr)
Create a Form from UFC input.
Definition: utils.h:228
IntegralType
Type of integral.
Definition: Form.h:31
@ interior_facet
Interior facet.
@ exterior_facet
Exterior facet.
la::SparsityPattern create_sparsity_pattern(const mesh::Topology &topology, const std::array< std::reference_wrapper< const DofMap >, 2 > &dofmaps, const std::set< IntegralType > &integrals)
Create a sparsity pattern for a given form.
Definition: utils.cpp:30
ElementDofLayout create_element_dof_layout(const ufcx_dofmap &dofmap, const mesh::CellType cell_type, const std::vector< int > &parent_map={})
Create an ElementDofLayout from a ufcx_dofmap.
Definition: utils.cpp:73
FunctionSpace create_functionspace(const std::shared_ptr< mesh::Mesh > &mesh, const basix::FiniteElement &e, int bs, const std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> &reorder_fn=nullptr)
Create a FunctionSpace from a Basix element.
Definition: utils.cpp:205
DofMap create_dofmap(MPI_Comm comm, const ElementDofLayout &layout, mesh::Topology &topology, const std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> &reorder_fn, const FiniteElement &element)
Create a dof map on mesh.
Definition: utils.cpp:140
std::pair< std::vector< T >, int > allocate_coefficient_storage(const Form< T > &form, IntegralType integral_type, int id)
Allocate storage for coefficients of a pair (integral_type, id) from a fem::Form form.
Definition: utils.h:675
void pack_coefficients(const Form< T > &form, IntegralType integral_type, int id, const std::span< T > &c, int cstride)
Pack coefficients of a Form for a given integral type and domain id.
Definition: utils.h:738
Mesh data structures and algorithms on meshes.
Definition: DofMap.h:30
CellType
Cell type identifier.
Definition: cell_types.h:22