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