DOLFINx
DOLFINx C++ interface
Form.h
1// Copyright (C) 2019-2020 Garth N. Wells and Chris Richardson
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 "FunctionSpace.h"
10#include <algorithm>
11#include <array>
12#include <dolfinx/common/IndexMap.h>
13#include <dolfinx/mesh/Mesh.h>
14#include <dolfinx/mesh/MeshTags.h>
15#include <functional>
16#include <memory>
17#include <string>
18#include <tuple>
19#include <vector>
20
21namespace dolfinx::fem
22{
23
24template <typename T>
25class Constant;
26template <typename T>
27class Function;
28
30enum class IntegralType : std::int8_t
31{
32 cell = 0,
33 exterior_facet = 1,
34 interior_facet = 2,
35 vertex = 3
36};
37
61template <typename T>
62class Form
63{
64 template <typename X, typename = void>
65 struct scalar_value_type
66 {
67 typedef X value_type;
68 };
69 template <typename X>
70 struct scalar_value_type<X, std::void_t<typename X::value_type>>
71 {
72 typedef typename X::value_type value_type;
73 };
74 using scalar_value_type_t = typename scalar_value_type<T>::value_type;
75
76public:
93 Form(const std::vector<std::shared_ptr<const FunctionSpace>>& function_spaces,
94 const std::map<
96 std::pair<
97 std::vector<std::pair<
98 int, std::function<void(T*, const T*, const T*,
99 const scalar_value_type_t*,
100 const int*, const std::uint8_t*)>>>,
101 const mesh::MeshTags<int>*>>& integrals,
102 const std::vector<std::shared_ptr<const Function<T>>>& coefficients,
103 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
105 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
106 : _function_spaces(function_spaces), _coefficients(coefficients),
107 _constants(constants), _mesh(mesh),
108 _needs_facet_permutations(needs_facet_permutations)
109 {
110 // Extract _mesh from FunctionSpace, and check they are the same
111 if (!_mesh and !function_spaces.empty())
112 _mesh = function_spaces[0]->mesh();
113 for (const auto& V : function_spaces)
114 {
115 if (_mesh != V->mesh())
116 throw std::runtime_error("Incompatible mesh");
117 }
118 if (!_mesh)
119 throw std::runtime_error("No mesh could be associated with the Form.");
120
121 // Store kernels, looping over integrals by domain type (dimension)
122 for (auto& integral_type : integrals)
123 {
124 const IntegralType type = integral_type.first;
125 // Loop over integrals kernels and set domains
126 switch (type)
127 {
129 for (auto& integral : integral_type.second.first)
130 _cell_integrals.insert({integral.first, {integral.second, {}}});
131 break;
133 for (auto& integral : integral_type.second.first)
134 {
135 _exterior_facet_integrals.insert(
136 {integral.first, {integral.second, {}}});
137 }
138 break;
140 for (auto& integral : integral_type.second.first)
141 {
142 _interior_facet_integrals.insert(
143 {integral.first, {integral.second, {}}});
144 }
145 break;
146 }
147
148 if (integral_type.second.second)
149 {
150 assert(_mesh == integral_type.second.second->mesh());
151 set_domains(type, *integral_type.second.second);
152 }
153 }
154
155 // FIXME: do this neatly via a static function
156 // Set markers for default integrals
157 set_default_domains(*_mesh);
158 }
159
161 Form(const Form& form) = delete;
162
164 Form(Form&& form) = default;
165
167 virtual ~Form() = default;
168
172 int rank() const { return _function_spaces.size(); }
173
176 std::shared_ptr<const mesh::Mesh> mesh() const { return _mesh; }
177
180 const std::vector<std::shared_ptr<const FunctionSpace>>&
182 {
183 return _function_spaces;
184 }
185
191 const std::function<void(T*, const T*, const T*, const scalar_value_type_t*,
192 const int*, const std::uint8_t*)>&
193 kernel(IntegralType type, int i) const
194 {
195 switch (type)
196 {
198 return get_kernel_from_integrals(_cell_integrals, i);
200 return get_kernel_from_integrals(_exterior_facet_integrals, i);
202 return get_kernel_from_integrals(_interior_facet_integrals, i);
203 default:
204 throw std::runtime_error(
205 "Cannot access kernel. Integral type not supported.");
206 }
207 }
208
211 std::set<IntegralType> integral_types() const
212 {
213 std::set<IntegralType> set;
214 if (!_cell_integrals.empty())
215 set.insert(IntegralType::cell);
216 if (!_exterior_facet_integrals.empty())
218 if (!_interior_facet_integrals.empty())
220
221 return set;
222 }
223
228 {
229 switch (type)
230 {
232 return _cell_integrals.size();
234 return _exterior_facet_integrals.size();
236 return _interior_facet_integrals.size();
237 default:
238 throw std::runtime_error("Integral type not supported.");
239 }
240 }
241
248 std::vector<int> integral_ids(IntegralType type) const
249 {
250 std::vector<int> ids;
251 switch (type)
252 {
254 std::transform(_cell_integrals.cbegin(), _cell_integrals.cend(),
255 std::back_inserter(ids),
256 [](auto& integral) { return integral.first; });
257 break;
259 std::transform(_exterior_facet_integrals.cbegin(),
260 _exterior_facet_integrals.cend(), std::back_inserter(ids),
261 [](auto& integral) { return integral.first; });
262 break;
264 std::transform(_interior_facet_integrals.cbegin(),
265 _interior_facet_integrals.cend(), std::back_inserter(ids),
266 [](auto& integral) { return integral.first; });
267 break;
268 default:
269 throw std::runtime_error(
270 "Cannot return IDs. Integral type not supported.");
271 }
272
273 return ids;
274 }
275
280 const std::vector<std::int32_t>& cell_domains(int i) const
281 {
282 auto it = _cell_integrals.find(i);
283 if (it == _cell_integrals.end())
284 throw std::runtime_error("No mesh entities for requested domain index.");
285 return it->second.second;
286 }
287
293 const std::vector<std::int32_t>& exterior_facet_domains(int i) const
294 {
295 auto it = _exterior_facet_integrals.find(i);
296 if (it == _exterior_facet_integrals.end())
297 throw std::runtime_error("No mesh entities for requested domain index.");
298 return it->second.second;
299 }
300
308 const std::vector<std::int32_t>& interior_facet_domains(int i) const
309 {
310 auto it = _interior_facet_integrals.find(i);
311 if (it == _interior_facet_integrals.end())
312 throw std::runtime_error("No mesh entities for requested domain index.");
313 return it->second.second;
314 }
315
317 const std::vector<std::shared_ptr<const Function<T>>>& coefficients() const
318 {
319 return _coefficients;
320 }
321
325 bool needs_facet_permutations() const { return _needs_facet_permutations; }
326
330 std::vector<int> coefficient_offsets() const
331 {
332 std::vector<int> n = {0};
333 for (const auto& c : _coefficients)
334 {
335 if (!c)
336 throw std::runtime_error("Not all form coefficients have been set.");
337 n.push_back(n.back() + c->function_space()->element()->space_dimension());
338 }
339 return n;
340 }
341
343 const std::vector<std::shared_ptr<const Constant<T>>>& constants() const
344 {
345 return _constants;
346 }
347
349 using scalar_type = T;
350
351private:
352 using kern
353 = std::function<void(T*, const T*, const T*, const scalar_value_type_t*,
354 const int*, const std::uint8_t*)>;
355
356 // Helper function to get the kernel for integral i from a map
357 // of integrals i.e. from _cell_integrals
358 // @param[in] integrals Map of integrals
359 // @param[in] i Domain index
360 // @return Function to call for tabulate_tensor
361 template <typename U>
362 const std::function<void(T*, const T*, const T*, const scalar_value_type_t*,
363 const int*, const std::uint8_t*)>&
364 get_kernel_from_integrals(const U& integrals, int i) const
365 {
366 auto it = integrals.find(i);
367 if (it == integrals.end())
368 throw std::runtime_error("No kernel for requested domain index.");
369 return it->second.first;
370 }
371
372 // Helper function to get a std::vector of (cell, local_facet) pairs
373 // corresponding to a given facet index.
374 // @param[in] f Facet index
375 // @param[in] f_to_c Facet to cell connectivity
376 // @param[in] c_to_f Cell to facet connectivity
377 // @return Vector of (cell, local_facet) pairs
378 template <int num_cells>
379 static std::array<std::array<std::int32_t, 2>, num_cells>
380 get_cell_local_facet_pairs(
381 std::int32_t f, const std::span<const std::int32_t>& cells,
383 {
384 // Loop over cells sharing facet
385 assert(cells.size() == num_cells);
386 std::array<std::array<std::int32_t, 2>, num_cells> cell_local_facet_pairs;
387 for (int c = 0; c < num_cells; ++c)
388 {
389 // Get local index of facet with respect to the cell
390 std::int32_t cell = cells[c];
391 auto cell_facets = c_to_f.links(cell);
392 auto facet_it = std::find(cell_facets.begin(), cell_facets.end(), f);
393 assert(facet_it != cell_facets.end());
394 int local_f = std::distance(cell_facets.begin(), facet_it);
395 cell_local_facet_pairs[c] = {cell, local_f};
396 }
397
398 return cell_local_facet_pairs;
399 }
400
401 // Set cell domains
402 template <typename iterator>
403 void set_cell_domains(
404 std::map<int, std::pair<kern, std::vector<std::int32_t>>>& integrals,
405 const iterator& tagged_cells_begin, const iterator& tagged_cells_end,
406 const std::vector<int>& tags)
407 {
408 // For cell integrals use all markers (but not on ghost entities)
409 for (auto c = tagged_cells_begin; c != tagged_cells_end; ++c)
410 {
411 const std::size_t pos = std::distance(tagged_cells_begin, c);
412 if (auto it = integrals.find(tags[pos]); it != integrals.end())
413 it->second.second.push_back(*c);
414 }
415 }
416
417 // Set exterior facet domains
418 template <typename iterator>
419 void set_exterior_facet_domains(
420 const mesh::Topology& topology,
421 std::map<int, std::pair<kern, std::vector<std::int32_t>>>& integrals,
422 const iterator& tagged_facets_begin, const iterator& tagged_facets_end,
423 const std::vector<int>& tags)
424 {
425 // When a mesh is not ghosted by cell, it is not straightforward
426 // to distinguish between (i) exterior facets and (ii) interior
427 // facets that are on a partition boundary. If there are no
428 // ghost cells, build a set of owned facts that are ghosted on
429 // another process to help determine if a facet is on an
430 // exterior boundary.
431 int tdim = topology.dim();
432 assert(topology.index_map(tdim));
433 assert(topology.index_map(tdim - 1));
434 const std::vector<std::int32_t> fwd_shared_facets
435 = topology.index_map(tdim)->overlapped()
436 ? std::vector<std::int32_t>()
437 : topology.index_map(tdim - 1)->shared_indices();
438
439 auto f_to_c = topology.connectivity(tdim - 1, tdim);
440 assert(f_to_c);
441 auto c_to_f = topology.connectivity(tdim, tdim - 1);
442 assert(c_to_f);
443 for (auto f = tagged_facets_begin; f != tagged_facets_end; ++f)
444 {
445 // All "owned" facets connected to one cell, that are not
446 // shared, should be external
447 // TODO: Consider removing this check and integrating over all
448 // tagged facets. This may be useful in a few cases.
449 if (f_to_c->num_links(*f) == 1)
450 {
451 if (!std::binary_search(fwd_shared_facets.begin(),
452 fwd_shared_facets.end(), *f))
453 {
454 const std::size_t pos = std::distance(tagged_facets_begin, f);
455 if (auto it = integrals.find(tags[pos]); it != integrals.end())
456 {
457 // There will only be one pair for an exterior facet integral
458 const std::array<std::int32_t, 2> pair
459 = get_cell_local_facet_pairs<1>(*f, f_to_c->links(*f),
460 *c_to_f)[0];
461 it->second.second.insert(it->second.second.end(), pair.cbegin(),
462 pair.cend());
463 }
464 }
465 }
466 }
467 }
468
469 // Set interior facet domains
470 template <typename iterator>
471 static void set_interior_facet_domains(
472 const mesh::Topology& topology,
473 std::map<int, std::pair<kern, std::vector<std::int32_t>>>& integrals,
474 const iterator& tagged_facets_begin, const iterator& tagged_facets_end,
475 const std::vector<int>& tags)
476 {
477 int tdim = topology.dim();
478 auto f_to_c = topology.connectivity(tdim - 1, tdim);
479 assert(f_to_c);
480 auto c_to_f = topology.connectivity(tdim, tdim - 1);
481 assert(c_to_f);
482 for (auto f = tagged_facets_begin; f != tagged_facets_end; ++f)
483 {
484 if (f_to_c->num_links(*f) == 2)
485 {
486 const std::size_t pos = std::distance(tagged_facets_begin, f);
487 if (auto it = integrals.find(tags[pos]); it != integrals.end())
488 {
489 const std::array<std::array<std::int32_t, 2>, 2> pairs
490 = get_cell_local_facet_pairs<2>(*f, f_to_c->links(*f), *c_to_f);
491 it->second.second.insert(it->second.second.end(), pairs[0].cbegin(),
492 pairs[0].cend());
493 it->second.second.insert(it->second.second.end(), pairs[1].cbegin(),
494 pairs[1].cend());
495 }
496 }
497 }
498 }
499
500 // Sets the entity indices to assemble over for kernels with a domain
501 // ID
502 // @param[in] type Integral type
503 // @param[in] marker MeshTags with domain ID. Entities with marker 'i'
504 // will be assembled over using the kernel with ID 'i'. The MeshTags
505 // is not stored.
506 void set_domains(IntegralType type, const mesh::MeshTags<int>& marker)
507 {
508 std::shared_ptr<const mesh::Mesh> mesh = marker.mesh();
509 const mesh::Topology& topology = mesh->topology();
510 const int tdim = topology.dim();
511 int dim = type == IntegralType::cell ? tdim : tdim - 1;
512 if (dim != marker.dim())
513 {
514 throw std::runtime_error("Invalid MeshTags dimension: "
515 + std::to_string(marker.dim()));
516 }
517
518 // Get mesh tag data
519 const std::vector<int>& tags = marker.values();
520 const std::vector<std::int32_t>& tagged_entities = marker.indices();
521 assert(topology.index_map(dim));
522 const auto entity_end
523 = std::lower_bound(tagged_entities.begin(), tagged_entities.end(),
524 topology.index_map(dim)->size_local());
525 switch (type)
526 {
528 set_cell_domains(_cell_integrals, tagged_entities.cbegin(), entity_end,
529 tags);
530 break;
531 default:
532 mesh->topology_mutable().create_connectivity(dim, tdim);
533 mesh->topology_mutable().create_connectivity(tdim, dim);
534 switch (type)
535 {
537 set_exterior_facet_domains(topology, _exterior_facet_integrals,
538 tagged_entities.cbegin(), entity_end, tags);
539 break;
541 set_interior_facet_domains(topology, _interior_facet_integrals,
542 tagged_entities.cbegin(), entity_end, tags);
543 break;
544 default:
545 throw std::runtime_error(
546 "Cannot set domains. Integral type not supported.");
547 }
548 }
549 }
550
556 void set_default_domains(const mesh::Mesh& mesh)
557 {
558 const mesh::Topology& topology = mesh.topology();
559 const int tdim = topology.dim();
560
561 // Cells. If there is a default integral, define it on all owned
562 // cells
563 for (auto& [domain_id, kernel_cells] : _cell_integrals)
564 {
565 if (domain_id == -1)
566 {
567 std::vector<std::int32_t>& cells = kernel_cells.second;
568 const int num_cells = topology.index_map(tdim)->size_local();
569 cells.resize(num_cells);
570 std::iota(cells.begin(), cells.end(), 0);
571 }
572 }
573
574 // Exterior facets. If there is a default integral, define it only
575 // on owned surface facets.
576
577 if (!_exterior_facet_integrals.empty())
578 {
579 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
580 mesh.topology_mutable().create_connectivity(tdim, tdim - 1);
581 }
582 const std::vector<std::int32_t> boundary_facets
583 = _exterior_facet_integrals.empty()
584 ? std::vector<std::int32_t>()
585 : mesh::exterior_facet_indices(topology);
586 for (auto& [domain_id, kernel_facets] : _exterior_facet_integrals)
587 {
588 if (domain_id == -1)
589 {
590 std::vector<std::int32_t>& facets = kernel_facets.second;
591 facets.clear();
592
593 auto f_to_c = topology.connectivity(tdim - 1, tdim);
594 assert(f_to_c);
595 auto c_to_f = topology.connectivity(tdim, tdim - 1);
596 assert(c_to_f);
597 for (std::int32_t f : boundary_facets)
598 {
599 // There will only be one pair for an exterior facet integral
600 std::array<std::int32_t, 2> pair
601 = get_cell_local_facet_pairs<1>(f, f_to_c->links(f), *c_to_f)[0];
602 facets.insert(facets.end(), pair.cbegin(), pair.cend());
603 }
604 }
605 }
606
607 // Interior facets. If there is a default integral, define it only on
608 // owned interior facets.
609 for (auto& [domain_id, kernel_facets] : _interior_facet_integrals)
610 {
611 if (domain_id == -1)
612 {
613 std::vector<std::int32_t>& facets = kernel_facets.second;
614 facets.clear();
615
616 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
617 auto f_to_c = topology.connectivity(tdim - 1, tdim);
618 assert(f_to_c);
619 mesh.topology_mutable().create_connectivity(tdim, tdim - 1);
620 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
621 assert(c_to_f);
622
623 // Get number of facets owned by this process
624 assert(topology.index_map(tdim - 1));
625 const int num_facets = topology.index_map(tdim - 1)->size_local();
626 facets.reserve(num_facets);
627 for (int f = 0; f < num_facets; ++f)
628 {
629 if (f_to_c->num_links(f) == 2)
630 {
631 const std::array<std::array<std::int32_t, 2>, 2> pairs
632 = get_cell_local_facet_pairs<2>(f, f_to_c->links(f), *c_to_f);
633 facets.insert(facets.end(), pairs[0].cbegin(), pairs[0].cend());
634 facets.insert(facets.end(), pairs[1].cbegin(), pairs[1].cend());
635 }
636 }
637 }
638 }
639 }
640
641 // Function spaces (one for each argument)
642 std::vector<std::shared_ptr<const FunctionSpace>> _function_spaces;
643
644 // Form coefficients
645 std::vector<std::shared_ptr<const Function<T>>> _coefficients;
646
647 // Constants associated with the Form
648 std::vector<std::shared_ptr<const Constant<T>>> _constants;
649
650 // The mesh
651 std::shared_ptr<const mesh::Mesh> _mesh;
652
653 // Cell integrals
654 std::map<int, std::pair<kern, std::vector<std::int32_t>>> _cell_integrals;
655
656 // Exterior facet integrals
657 std::map<int, std::pair<kern, std::vector<std::int32_t>>>
658 _exterior_facet_integrals;
659
660 // Interior facet integrals
661 std::map<int, std::pair<kern, std::vector<std::int32_t>>>
662 _interior_facet_integrals;
663
664 // True if permutation data needs to be passed into these integrals
665 bool _needs_facet_permutations;
666};
667} // namespace dolfinx::fem
Constant value which can be attached to a Form. Constants may be scalar (rank 0), vector (rank 1),...
Definition: Constant.h:20
A representation of finite element variational forms.
Definition: Form.h:63
bool needs_facet_permutations() const
Get bool indicating whether permutation data needs to be passed into these integrals.
Definition: Form.h:325
Form(Form &&form)=default
Move constructor.
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
Form(const Form &form)=delete
Copy constructor.
virtual ~Form()=default
Destructor.
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::function< void(T *, const T *, const T *, const scalar_value_type_t *, const int *, const std::uint8_t *)> & kernel(IntegralType type, int i) const
Get the function for 'kernel' for integral i of given type.
Definition: Form.h:193
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
const std::vector< std::shared_ptr< const Constant< T > > > & constants() const
Access constants.
Definition: Form.h:343
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
Form(const std::vector< std::shared_ptr< const FunctionSpace > > &function_spaces, const std::map< IntegralType, std::pair< std::vector< std::pair< int, std::function< void(T *, const T *, const T *, const scalar_value_type_t *, const int *, const std::uint8_t *)> > >, const mesh::MeshTags< int > * > > &integrals, const std::vector< std::shared_ptr< const Function< T > > > &coefficients, const std::vector< std::shared_ptr< const Constant< T > > > &constants, bool needs_facet_permutations, const std::shared_ptr< const mesh::Mesh > &mesh=nullptr)
Create a finite element form.
Definition: Form.h:93
int num_integrals(IntegralType type) const
Number of integrals of given type.
Definition: Form.h:227
T scalar_type
Scalar type (T)
Definition: Form.h:349
This class represents a function in a finite element function space , given by.
Definition: Function.h:45
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:26
std::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:111
MeshTags associate values with mesh entities.
Definition: MeshTags.h:36
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
Finite element method functionality.
Definition: assemble_matrix_impl.h:25
IntegralType
Type of integral.
Definition: Form.h:31
@ interior_facet
Interior facet.
@ exterior_facet
Exterior facet.
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:570