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 <dolfinx/la/SparsityPattern.h>
17#include <dolfinx/mesh/cell_types.h>
18#include <functional>
19#include <memory>
20#include <set>
21#include <stdexcept>
22#include <string>
23#include <ufcx.h>
24#include <utility>
25#include <vector>
26#include <xtensor/xadapt.hpp>
27#include <xtensor/xtensor.hpp>
28#include <xtl/xspan.hpp>
29
32
33namespace basix
34{
35class FiniteElement;
36}
37
38namespace dolfinx::common
39{
40class IndexMap;
41}
42
43namespace dolfinx::mesh
44{
45class Mesh;
46class Topology;
47} // namespace dolfinx::mesh
48
49namespace dolfinx::fem
50{
51class FunctionSpace;
52
60template <typename T>
61std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace>, 2>>>
62extract_function_spaces(const std::vector<std::vector<const Form<T>*>>& a)
63{
64 std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace>, 2>>>
65 spaces(a.size(),
66 std::vector<std::array<std::shared_ptr<const FunctionSpace>, 2>>(
67 a[0].size()));
68 for (std::size_t i = 0; i < a.size(); ++i)
69 {
70 for (std::size_t j = 0; j < a[i].size(); ++j)
71 {
72 if (const Form<T>* form = a[i][j]; form)
73 spaces[i][j] = {form->function_spaces()[0], form->function_spaces()[1]};
74 }
75 }
76 return spaces;
77}
78
83 const mesh::Topology& topology,
84 const std::array<std::reference_wrapper<const DofMap>, 2>& dofmaps,
85 const std::set<IntegralType>& integrals);
86
92template <typename T>
94{
95 if (a.rank() != 2)
96 {
97 throw std::runtime_error(
98 "Cannot create sparsity pattern. Form is not a bilinear form");
99 }
100
101 // Get dof maps and mesh
102 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
103 *a.function_spaces().at(0)->dofmap(),
104 *a.function_spaces().at(1)->dofmap()};
105 std::shared_ptr mesh = a.mesh();
106 assert(mesh);
107
108 const std::set<IntegralType> types = a.integral_types();
109 if (types.find(IntegralType::interior_facet) != types.end()
110 or types.find(IntegralType::exterior_facet) != types.end())
111 {
112 // FIXME: cleanup these calls? Some of the happen internally again.
113 const int tdim = mesh->topology().dim();
114 mesh->topology_mutable().create_entities(tdim - 1);
115 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
116 }
117
118 return create_sparsity_pattern(mesh->topology(), dofmaps, types);
119}
120
122ElementDofLayout create_element_dof_layout(const ufcx_dofmap& dofmap,
123 const mesh::CellType cell_type,
124 const std::vector<int>& parent_map
125 = {});
126
135DofMap
136create_dofmap(MPI_Comm comm, const ElementDofLayout& layout,
137 mesh::Topology& topology,
138 const std::function<std::vector<int>(
139 const graph::AdjacencyList<std::int32_t>&)>& reorder_fn,
140 const FiniteElement& element);
141
145std::vector<std::string> get_coefficient_names(const ufcx_form& ufcx_form);
146
150std::vector<std::string> get_constant_names(const ufcx_form& ufcx_form);
151
159template <typename T>
161 const ufcx_form& ufcx_form,
162 const std::vector<std::shared_ptr<const FunctionSpace>>& spaces,
163 const std::vector<std::shared_ptr<const Function<T>>>& coefficients,
164 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
165 const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
166 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
167{
168 if (ufcx_form.rank != (int)spaces.size())
169 throw std::runtime_error("Wrong number of argument spaces for Form.");
170 if (ufcx_form.num_coefficients != (int)coefficients.size())
171 {
172 throw std::runtime_error(
173 "Mismatch between number of expected and provided Form coefficients.");
174 }
175 if (ufcx_form.num_constants != (int)constants.size())
176 {
177 throw std::runtime_error(
178 "Mismatch between number of expected and provided Form constants.");
179 }
180
181 // Check argument function spaces
182#ifndef NDEBUG
183 for (std::size_t i = 0; i < spaces.size(); ++i)
184 {
185 assert(spaces[i]->element());
186 ufcx_finite_element* ufcx_element = ufcx_form.finite_elements[i];
187 assert(ufcx_element);
188 if (std::string(ufcx_element->signature)
189 != spaces[i]->element()->signature())
190 {
191 throw std::runtime_error(
192 "Cannot create form. Wrong type of function space for argument.");
193 }
194 }
195#endif
196
197 // Get list of integral IDs, and load tabulate tensor into memory for
198 // each
199 using kern = std::function<void(T*, const T*, const T*, const double*,
200 const int*, const std::uint8_t*)>;
201 std::map<IntegralType, std::pair<std::vector<std::pair<int, kern>>,
202 const mesh::MeshTags<int>*>>
203 integral_data;
204
205 bool needs_facet_permutations = false;
206
207 // Attach cell kernels
208 std::vector<int> cell_integral_ids(ufcx_form.integral_ids(cell),
209 ufcx_form.integral_ids(cell)
210 + ufcx_form.num_integrals(cell));
211 for (int i = 0; i < ufcx_form.num_integrals(cell); ++i)
212 {
213 ufcx_integral* integral = ufcx_form.integrals(cell)[i];
214 assert(integral);
215
216 kern k = nullptr;
217 if constexpr (std::is_same<T, float>::value)
218 k = integral->tabulate_tensor_float32;
219 else if constexpr (std::is_same<T, std::complex<float>>::value)
220 {
221 k = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
222 const int*, const unsigned char*)>(
223 integral->tabulate_tensor_complex64);
224 }
225 else if constexpr (std::is_same<T, double>::value)
226 k = integral->tabulate_tensor_float64;
227 else if constexpr (std::is_same<T, std::complex<double>>::value)
228 {
229 k = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
230 const int*, const unsigned char*)>(
231 integral->tabulate_tensor_complex128);
232 }
233 assert(k);
234
235 integral_data[IntegralType::cell].first.emplace_back(cell_integral_ids[i],
236 k);
237 if (integral->needs_facet_permutations)
238 needs_facet_permutations = true;
239 }
240
241 // Attach cell subdomain data
242 if (auto it = subdomains.find(IntegralType::cell);
243 it != subdomains.end() and !cell_integral_ids.empty())
244 {
245 integral_data[IntegralType::cell].second = it->second;
246 }
247
248 // FIXME: Can facets be handled better?
249
250 // Create facets, if required
251 if (ufcx_form.num_integrals(exterior_facet) > 0
252 or ufcx_form.num_integrals(interior_facet) > 0)
253 {
254 if (!spaces.empty())
255 {
256 auto mesh = spaces[0]->mesh();
257 const int tdim = mesh->topology().dim();
258 spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
259 }
260 }
261
262 // Attach exterior facet kernels
263 std::vector<int> exterior_facet_integral_ids(
264 ufcx_form.integral_ids(exterior_facet),
265 ufcx_form.integral_ids(exterior_facet)
266 + ufcx_form.num_integrals(exterior_facet));
267 for (int i = 0; i < ufcx_form.num_integrals(exterior_facet); ++i)
268 {
269 ufcx_integral* integral = ufcx_form.integrals(exterior_facet)[i];
270 assert(integral);
271
272 kern k = nullptr;
273 if constexpr (std::is_same<T, float>::value)
274 k = integral->tabulate_tensor_float32;
275 else if constexpr (std::is_same<T, std::complex<float>>::value)
276 {
277 k = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
278 const int*, const unsigned char*)>(
279 integral->tabulate_tensor_complex64);
280 }
281 else if constexpr (std::is_same<T, double>::value)
282 k = integral->tabulate_tensor_float64;
283 else if constexpr (std::is_same<T, std::complex<double>>::value)
284 {
285 k = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
286 const int*, const unsigned char*)>(
287 integral->tabulate_tensor_complex128);
288 }
289 assert(k);
290
291 integral_data[IntegralType::exterior_facet].first.emplace_back(
292 exterior_facet_integral_ids[i], k);
293 if (integral->needs_facet_permutations)
294 needs_facet_permutations = true;
295 }
296
297 // Attach exterior facet subdomain data
298 if (auto it = subdomains.find(IntegralType::exterior_facet);
299 it != subdomains.end() and !exterior_facet_integral_ids.empty())
300 {
301 integral_data[IntegralType::exterior_facet].second = it->second;
302 }
303
304 // Attach interior facet kernels
305 std::vector<int> interior_facet_integral_ids(
306 ufcx_form.integral_ids(interior_facet),
307 ufcx_form.integral_ids(interior_facet)
308 + ufcx_form.num_integrals(interior_facet));
309 for (int i = 0; i < ufcx_form.num_integrals(interior_facet); ++i)
310 {
311 ufcx_integral* integral = ufcx_form.integrals(interior_facet)[i];
312 assert(integral);
313
314 kern k = nullptr;
315 if constexpr (std::is_same<T, float>::value)
316 k = integral->tabulate_tensor_float32;
317 else if constexpr (std::is_same<T, std::complex<float>>::value)
318 {
319 k = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
320 const int*, const unsigned char*)>(
321 integral->tabulate_tensor_complex64);
322 }
323 else if constexpr (std::is_same<T, double>::value)
324 k = integral->tabulate_tensor_float64;
325 else if constexpr (std::is_same<T, std::complex<double>>::value)
326 {
327 k = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
328 const int*, const unsigned char*)>(
329 integral->tabulate_tensor_complex128);
330 }
331 assert(k);
332
333 integral_data[IntegralType::interior_facet].first.emplace_back(
334 interior_facet_integral_ids[i], k);
335 if (integral->needs_facet_permutations)
336 needs_facet_permutations = true;
337 }
338
339 // Attach interior facet subdomain data
340 if (auto it = subdomains.find(IntegralType::interior_facet);
341 it != subdomains.end() and !interior_facet_integral_ids.empty())
342 {
343 integral_data[IntegralType::interior_facet].second = it->second;
344 }
345
346 return Form<T>(spaces, integral_data, coefficients, constants,
347 needs_facet_permutations, mesh);
348}
349
359template <typename T>
361 const ufcx_form& ufcx_form,
362 const std::vector<std::shared_ptr<const FunctionSpace>>& spaces,
363 const std::map<std::string, std::shared_ptr<const Function<T>>>&
364 coefficients,
365 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
366 const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
367 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
368{
369 // Place coefficients in appropriate order
370 std::vector<std::shared_ptr<const Function<T>>> coeff_map;
371 for (const std::string& name : get_coefficient_names(ufcx_form))
372 {
373 if (auto it = coefficients.find(name); it != coefficients.end())
374 coeff_map.push_back(it->second);
375 else
376 {
377 throw std::runtime_error("Form coefficient \"" + name
378 + "\" not provided.");
379 }
380 }
381
382 // Place constants in appropriate order
383 std::vector<std::shared_ptr<const Constant<T>>> const_map;
384 for (const std::string& name : get_constant_names(ufcx_form))
385 {
386 if (auto it = constants.find(name); it != constants.end())
387 const_map.push_back(it->second);
388 else
389 throw std::runtime_error("Form constant \"" + name + "\" not provided.");
390 }
391
392 return create_form(ufcx_form, spaces, coeff_map, const_map, subdomains, mesh);
393}
394
406template <typename T>
408 ufcx_form* (*fptr)(),
409 const std::vector<std::shared_ptr<const FunctionSpace>>& spaces,
410 const std::map<std::string, std::shared_ptr<const Function<T>>>&
411 coefficients,
412 const std::map<std::string, std::shared_ptr<const Constant<T>>>& constants,
413 const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
414 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
415{
416 ufcx_form* form = fptr();
417 Form<T> L = create_form<T>(*form, spaces, coefficients, constants, subdomains,
418 mesh);
419 std::free(form);
420 return L;
421}
422
431FunctionSpace
432create_functionspace(const std::shared_ptr<mesh::Mesh>& mesh,
433 const basix::FiniteElement& e, int bs,
434 const std::function<std::vector<int>(
435 const graph::AdjacencyList<std::int32_t>&)>& reorder_fn
436 = nullptr);
437
450FunctionSpace create_functionspace(
451 ufcx_function_space* (*fptr)(const char*), const std::string& function_name,
452 const std::shared_ptr<mesh::Mesh>& mesh,
453 const std::function<
454 std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>& reorder_fn
455 = nullptr);
456
458namespace impl
459{
461template <typename T>
462xtl::span<const std::uint32_t> get_cell_orientation_info(
463 const std::vector<std::shared_ptr<const Function<T>>>& coefficients)
464{
465 bool needs_dof_transformations = false;
466 for (auto coeff : coefficients)
467 {
468 std::shared_ptr<const FiniteElement> element
469 = coeff->function_space()->element();
470 if (element->needs_dof_transformations())
471 {
472 needs_dof_transformations = true;
473 break;
474 }
475 }
476
477 xtl::span<const std::uint32_t> cell_info;
478 if (needs_dof_transformations)
479 {
480 auto mesh = coefficients.front()->function_space()->mesh();
481 mesh->topology_mutable().create_entity_permutations();
482 cell_info = xtl::span(mesh->topology().get_cell_permutation_info());
483 }
484
485 return cell_info;
486}
487
488// Pack a single coefficient for a single cell
489template <typename T, int _bs, typename Functor>
490static inline void pack(const xtl::span<T>& coeffs, std::int32_t cell, int bs,
491 const xtl::span<const T>& v,
492 const xtl::span<const std::uint32_t>& cell_info,
493 const DofMap& dofmap, Functor transform)
494{
495 auto dofs = dofmap.cell_dofs(cell);
496 for (std::size_t i = 0; i < dofs.size(); ++i)
497 {
498 if constexpr (_bs < 0)
499 {
500 const int pos_c = bs * i;
501 const int pos_v = bs * dofs[i];
502 for (int k = 0; k < bs; ++k)
503 coeffs[pos_c + k] = v[pos_v + k];
504 }
505 else
506 {
507 const int pos_c = _bs * i;
508 const int pos_v = _bs * dofs[i];
509 for (int k = 0; k < _bs; ++k)
510 coeffs[pos_c + k] = v[pos_v + k];
511 }
512 }
513
514 transform(coeffs, cell_info, cell, 1);
515}
516
530template <typename T, typename E, typename Functor>
531void pack_coefficient_entity(const xtl::span<T>& c, int cstride,
532 const Function<T>& u,
533 const xtl::span<const std::uint32_t>& cell_info,
534 const E& entities, Functor fetch_cells,
535 std::int32_t offset)
536{
537 // Read data from coefficient "u"
538 const xtl::span<const T>& v = u.x()->array();
539 const DofMap& dofmap = *u.function_space()->dofmap();
540 std::shared_ptr<const FiniteElement> element = u.function_space()->element();
541 int space_dim = element->space_dimension();
542 const auto transformation
543 = element->get_dof_transformation_function<T>(false, true);
544
545 const int bs = dofmap.bs();
546 switch (bs)
547 {
548 case 1:
549 for (std::size_t e = 0; e < entities.size(); ++e)
550 {
551 std::int32_t cell = fetch_cells(entities[e]);
552 auto cell_coeff = c.subspan(e * cstride + offset, space_dim);
553 pack<T, 1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
554 }
555 break;
556 case 2:
557 for (std::size_t e = 0; e < entities.size(); ++e)
558 {
559 std::int32_t cell = fetch_cells(entities[e]);
560 auto cell_coeff = c.subspan(e * cstride + offset, space_dim);
561 pack<T, 2>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
562 }
563 break;
564 case 3:
565 for (std::size_t e = 0; e < entities.size(); ++e)
566 {
567 std::int32_t cell = fetch_cells(entities[e]);
568 auto cell_coeff = c.subspan(e * cstride + offset, space_dim);
569 pack<T, 3>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
570 }
571 break;
572 default:
573 for (std::size_t e = 0; e < entities.size(); ++e)
574 {
575 std::int32_t cell = fetch_cells(entities[e]);
576 auto cell_coeff = c.subspan(e * cstride + offset, space_dim);
577 pack<T, -1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation);
578 }
579 break;
580 }
581}
582
583} // namespace impl
584
591template <typename T>
592std::pair<std::vector<T>, int>
594 int id)
595{
596 // Get form coefficient offsets and dofmaps
597 const std::vector<std::shared_ptr<const Function<T>>>& coefficients
598 = form.coefficients();
599 const std::vector<int> offsets = form.coefficient_offsets();
600
601 std::size_t num_entities = 0;
602 int cstride = 0;
603 if (!coefficients.empty())
604 {
605 cstride = offsets.back();
606 switch (integral_type)
607 {
609 num_entities = form.cell_domains(id).size();
610 break;
612 num_entities = form.exterior_facet_domains(id).size();
613 break;
615 num_entities = form.interior_facet_domains(id).size() * 2;
616 break;
617 default:
618 throw std::runtime_error(
619 "Could not allocate coefficient data. Integral type not supported.");
620 }
621 }
622
623 return {std::vector<T>(num_entities * cstride), cstride};
624}
625
630template <typename T>
631std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>
633{
634 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>> coeffs;
635 for (auto integral_type : form.integral_types())
636 {
637 for (int id : form.integral_ids(integral_type))
638 {
639 coeffs.emplace_hint(
640 coeffs.end(), std::pair(integral_type, id),
641 allocate_coefficient_storage(form, integral_type, id));
642 }
643 }
644
645 return coeffs;
646}
647
655template <typename T>
656void pack_coefficients(const Form<T>& form, IntegralType integral_type, int id,
657 const xtl::span<T>& c, int cstride)
658{
659 // Get form coefficient offsets and dofmaps
660 const std::vector<std::shared_ptr<const Function<T>>>& coefficients
661 = form.coefficients();
662 const std::vector<int> offsets = form.coefficient_offsets();
663
664 if (!coefficients.empty())
665 {
666 xtl::span<const std::uint32_t> cell_info
667 = impl::get_cell_orientation_info(coefficients);
668
669 switch (integral_type)
670 {
672 {
673 auto fetch_cell = [](auto entity) { return entity; };
674 const std::vector<std::int32_t>& cells = form.cell_domains(id);
675 // Iterate over coefficients
676 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
677 {
678 impl::pack_coefficient_entity(c, cstride, *coefficients[coeff],
679 cell_info, cells, fetch_cell,
680 offsets[coeff]);
681 }
682 break;
683 }
685 {
686 const std::vector<std::pair<std::int32_t, int>>& facets
687 = form.exterior_facet_domains(id);
688
689 // Create lambda function fetching cell index from exterior facet entity
690 auto fetch_cell = [](auto& entity) { return entity.first; };
691
692 // Iterate over coefficients
693 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
694 {
695 impl::pack_coefficient_entity(c, cstride, *coefficients[coeff],
696 cell_info, facets, fetch_cell,
697 offsets[coeff]);
698 }
699
700 break;
701 }
703 {
704 const std::vector<std::tuple<std::int32_t, int, std::int32_t, int>>&
705 facets
706 = form.interior_facet_domains(id);
707 // Lambda functions to fetch cell index from interior facet entity
708 auto fetch_cell0 = [](auto& entity) { return std::get<0>(entity); };
709 auto fetch_cell1 = [](auto& entity) { return std::get<2>(entity); };
710
711 // Iterate over coefficients
712 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
713 {
714 // Pack coefficient ['+']
715 impl::pack_coefficient_entity(c, 2 * cstride, *coefficients[coeff],
716 cell_info, facets, fetch_cell0,
717 2 * offsets[coeff]);
718 // Pack coefficient ['-']
719 impl::pack_coefficient_entity(c, 2 * cstride, *coefficients[coeff],
720 cell_info, facets, fetch_cell1,
721 offsets[coeff] + offsets[coeff + 1]);
722 }
723 break;
724 }
725 default:
726 throw std::runtime_error(
727 "Could not pack coefficient. Integral type not supported.");
728 }
729 }
730}
731
733template <typename T>
735 const ufcx_expression& expression,
736 const std::vector<std::shared_ptr<const fem::Function<T>>>& coefficients,
737 const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants,
738 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr,
739 const std::shared_ptr<const fem::FunctionSpace>& argument_function_space
740 = nullptr)
741{
742 if (expression.rank > 0 and !argument_function_space)
743 {
744 throw std::runtime_error(
745 "Expression has Argument but no Argument function space was provided.");
746 }
747
748 const int size = expression.num_points * expression.topological_dimension;
749 const xt::xtensor<double, 2> points = xt::adapt(
750 expression.points, size, xt::no_ownership(),
751 std::vector<std::size_t>(
752 {static_cast<std::size_t>(expression.num_points),
753 static_cast<std::size_t>(expression.topological_dimension)}));
754
755 std::vector<int> value_shape;
756 for (int i = 0; i < expression.num_components; ++i)
757 value_shape.push_back(expression.value_shape[i]);
758
759 std::function<void(T*, const T*, const T*, const double*, const int*,
760 const std::uint8_t*)>
761 tabulate_tensor = nullptr;
762 if constexpr (std::is_same<T, float>::value)
763 tabulate_tensor = expression.tabulate_tensor_float32;
764 else if constexpr (std::is_same<T, std::complex<float>>::value)
765 {
766 tabulate_tensor
767 = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
768 const int*, const unsigned char*)>(
769 expression.tabulate_tensor_complex64);
770 }
771 else if constexpr (std::is_same<T, double>::value)
772 tabulate_tensor = expression.tabulate_tensor_float64;
773 else if constexpr (std::is_same<T, std::complex<double>>::value)
774 {
775 tabulate_tensor
776 = reinterpret_cast<void (*)(T*, const T*, const T*, const double*,
777 const int*, const unsigned char*)>(
778 expression.tabulate_tensor_complex128);
779 }
780 else
781 {
782 throw std::runtime_error("Type not supported.");
783 }
784 assert(tabulate_tensor);
785
786 return fem::Expression(coefficients, constants, points, tabulate_tensor,
787 value_shape, mesh, argument_function_space);
788}
789
792template <typename T>
794 const ufcx_expression& expression,
795 const std::map<std::string, std::shared_ptr<const fem::Function<T>>>&
796 coefficients,
797 const std::map<std::string, std::shared_ptr<const fem::Constant<T>>>&
798 constants,
799 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr,
800 const std::shared_ptr<const fem::FunctionSpace>& argument_function_space
801 = nullptr)
802{
803 // Place coefficients in appropriate order
804 std::vector<std::shared_ptr<const fem::Function<T>>> coeff_map;
805 std::vector<std::string> coefficient_names;
806 for (int i = 0; i < expression.num_coefficients; ++i)
807 coefficient_names.push_back(expression.coefficient_names[i]);
808
809 for (const std::string& name : coefficient_names)
810 {
811 if (auto it = coefficients.find(name); it != coefficients.end())
812 coeff_map.push_back(it->second);
813 else
814 {
815 throw std::runtime_error("Expression coefficient \"" + name
816 + "\" not provided.");
817 }
818 }
819
820 // Place constants in appropriate order
821 std::vector<std::shared_ptr<const fem::Constant<T>>> const_map;
822 std::vector<std::string> constant_names;
823 for (int i = 0; i < expression.num_constants; ++i)
824 constant_names.push_back(expression.constant_names[i]);
825
826 for (const std::string& name : constant_names)
827 {
828 if (auto it = constants.find(name); it != constants.end())
829 const_map.push_back(it->second);
830 else
831 {
832 throw std::runtime_error("Expression constant \"" + name
833 + "\" not provided.");
834 }
835 }
836
837 return create_expression(expression, coeff_map, const_map, mesh,
838 argument_function_space);
839}
840
846template <typename T>
847void pack_coefficients(const Form<T>& form,
848 std::map<std::pair<IntegralType, int>,
849 std::pair<std::vector<T>, int>>& coeffs)
850{
851 for (auto& [key, val] : coeffs)
852 pack_coefficients<T>(form, key.first, key.second, val.first, val.second);
853}
854
861template <typename T>
862std::pair<std::vector<T>, int>
864 const xtl::span<const std::int32_t>& cells)
865{
866 // Get form coefficient offsets and dofmaps
867 const std::vector<std::shared_ptr<const Function<T>>>& coefficients
868 = u.coefficients();
869 const std::vector<int> offsets = u.coefficient_offsets();
870
871 // Copy data into coefficient array
872 const int cstride = offsets.back();
873 std::vector<T> c(cells.size() * offsets.back());
874 if (!coefficients.empty())
875 {
876 xtl::span<const std::uint32_t> cell_info
877 = impl::get_cell_orientation_info(coefficients);
878
879 // Iterate over coefficients
880 for (std::size_t coeff = 0; coeff < coefficients.size(); ++coeff)
881 impl::pack_coefficient_entity(
882 xtl::span(c), cstride, *coefficients[coeff], cell_info, cells,
883 [](std::int32_t entity) { return entity; }, offsets[coeff]);
884 }
885 return {std::move(c), cstride};
886}
887
890template <typename U>
891std::vector<typename U::scalar_type> pack_constants(const U& u)
892{
893 using T = typename U::scalar_type;
894 const std::vector<std::shared_ptr<const Constant<T>>>& constants
895 = u.constants();
896
897 // Calculate size of array needed to store packed constants
898 std::int32_t size = std::accumulate(constants.cbegin(), constants.cend(), 0,
899 [](std::int32_t sum, auto& constant)
900 { return sum + constant->value.size(); });
901
902 // Pack constants
903 std::vector<T> constant_values(size);
904 std::int32_t offset = 0;
905 for (auto& constant : constants)
906 {
907 const std::vector<T>& value = constant->value;
908 std::copy(value.cbegin(), value.cend(),
909 std::next(constant_values.begin(), offset));
910 offset += value.size();
911 }
912
913 return constant_values;
914}
915
916} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
Constant value which can be attached to a Form. Constants may be scalar (rank 0), vector (rank 1),...
Definition: Constant.h:21
Degree-of-freedom map.
Definition: DofMap.h:71
int bs() const noexcept
Return the block size for the dofmap.
Definition: DofMap.cpp:204
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:106
const std::vector< std::shared_ptr< const fem::Function< T > > > & coefficients() const
Get coefficients.
Definition: Expression.h:89
A representation of finite element variational forms.
Definition: Form.h:62
const std::vector< std::pair< std::int32_t, int > > & 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:281
int rank() const
Rank of the form (bilinear form = 2, linear form = 1, functional = 0, etc)
Definition: Form.h:160
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:268
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:319
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:236
std::set< IntegralType > integral_types() const
Get types of integrals in the form.
Definition: Form.h:199
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh for the form.
Definition: Form.h:164
const std::vector< std::shared_ptr< const fem::Function< T > > > & coefficients() const
Access coefficients.
Definition: Form.h:306
const std::vector< std::shared_ptr< const fem::FunctionSpace > > & function_spaces() const
Return function spaces for all arguments.
Definition: Form.h:169
const std::vector< std::tuple< std::int32_t, int, std::int32_t, int > > & interior_facet_domains(int i) const
Get the list of (cell_index_0, local_facet_index_0, cell_index_1, local_facet_index_1) tuples for the...
Definition: Form.h:296
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:46
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:57
void pack_coefficient_entity(const xtl::span< T > &c, int cstride, const Function< T > &u, const xtl::span< const std::uint32_t > &cell_info, const E &entities, Functor fetch_cells, std::int32_t offset)
Pack a single coefficient for a set of active entities.
Definition: utils.h:531
Miscellaneous classes, functions and types.
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
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:197
fem::Expression< T > create_expression(const ufcx_expression &expression, const std::vector< std::shared_ptr< const fem::Function< T > > > &coefficients, const std::vector< std::shared_ptr< const fem::Constant< T > > > &constants, const std::shared_ptr< const mesh::Mesh > &mesh=nullptr, const std::shared_ptr< const fem::FunctionSpace > &argument_function_space=nullptr)
Create Expression from UFC.
Definition: utils.h:734
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:188
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:62
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:891
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:160
IntegralType
Type of integral.
Definition: Form.h:30
@ 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:31
void pack_coefficients(const Form< T > &form, IntegralType integral_type, int id, const xtl::span< T > &c, int cstride)
Pack coefficients of a Form for a given integral type and domain id.
Definition: utils.h:656
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:74
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:206
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:141
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:593
Mesh data structures and algorithms on meshes.
Definition: DofMap.h:30
CellType
Cell type identifier.
Definition: cell_types.h:22