DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
Form.h
1// Copyright (C) 2019-2025 Garth N. Wells, Chris Richardson, Joseph P. Dean and
2// Jørgen S. Dokken
3//
4// This file is part of DOLFINx (https://www.fenicsproject.org)
5//
6// SPDX-License-Identifier: LGPL-3.0-or-later
7
8#pragma once
9
10#include "FunctionSpace.h"
11#include "traits.h"
12#include <algorithm>
13#include <basix/mdspan.hpp>
14#include <concepts>
15#include <cstdint>
16#include <dolfinx/common/IndexMap.h>
17#include <dolfinx/common/types.h>
18#include <dolfinx/mesh/EntityMap.h>
19#include <dolfinx/mesh/Mesh.h>
20#include <functional>
21#include <map>
22#include <memory>
23#include <optional>
24#include <ranges>
25#include <span>
26#include <tuple>
27#include <utility>
28#include <vector>
29
30namespace dolfinx::fem
31{
32template <dolfinx::scalar T>
33class Constant;
34template <dolfinx::scalar T, std::floating_point U>
35class Function;
36
38enum class IntegralType : std::int8_t
39{
40 cell = 0,
43 vertex = 3
44};
45
49template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
51{
57 template <typename K, typename V, typename W>
58 requires std::is_convertible_v<
59 std::remove_cvref_t<K>,
60 std::function<void(T*, const T*, const T*, const U*,
61 const int*, const uint8_t*, void*)>>
62 and std::is_convertible_v<std::remove_cvref_t<V>,
63 std::vector<std::int32_t>>
64 and std::is_convertible_v<std::remove_cvref_t<W>,
65 std::vector<int>>
67 : kernel(std::forward<K>(kernel)), entities(std::forward<V>(entities)),
68 coeffs(std::forward<W>(coeffs))
69 {
70 }
71
73 std::function<void(T*, const T*, const T*, const U*, const int*,
74 const uint8_t*, void*)>
76
79 std::vector<std::int32_t> entities;
80
83 std::vector<int> coeffs;
84};
85
114template <dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_t<T>>
115class Form
116{
117public:
119 using scalar_type = T;
120
122 using geometry_type = U;
123
154 template <typename X>
155 requires std::is_convertible_v<
156 std::remove_cvref_t<X>,
157 std::map<std::tuple<IntegralType, int, int>,
160 const std::vector<std::shared_ptr<const FunctionSpace<geometry_type>>>& V,
161 X&& integrals, std::shared_ptr<const mesh::Mesh<geometry_type>> mesh,
162 const std::vector<
163 std::shared_ptr<const Function<scalar_type, geometry_type>>>&
165 const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
166 constants,
168 const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
169 entity_maps)
170 : _function_spaces(V), _integrals(std::forward<X>(integrals)),
171 _mesh(mesh), _coefficients(coefficients), _constants(constants),
172 _needs_facet_permutations(needs_facet_permutations)
173 {
174 if (!_mesh)
175 throw std::runtime_error("Form Mesh is null.");
176
177 // A helper function to find the correct entity map for a given mesh
178 auto get_entity_map
179 = [mesh, &entity_maps](auto& mesh0) -> const mesh::EntityMap&
180 {
181 auto it = std::ranges::find_if(
182 entity_maps,
183 [mesh, mesh0](const mesh::EntityMap& em)
184 {
185 return ((em.topology() == mesh0->topology()
186 and em.sub_topology() == mesh->topology()))
187 or ((em.sub_topology() == mesh0->topology()
188 and em.topology() == mesh->topology()));
189 });
190
191 if (it == entity_maps.end())
192 {
193 throw std::runtime_error(
194 "Incompatible mesh. argument entity_maps must be provided.");
195 }
196 return *it;
197 };
198
199 // A helper function to compute the (cell, local_facet) pairs in the
200 // argument/coefficient domain from the (cell, local_facet) pairs in
201 // `this->mesh()`.
202 auto compute_facet_domains
203 = [&](const auto& int_ents_mesh, int codim, const auto& c_to_f,
204 const auto& emap, bool inverse)
205 {
206 // TODO: This function would be much neater using
207 // `std::views::stride(2)` from C++ 23
208
209 // Get a list of entities to map to the argument/coefficient
210 // domain
211 std::vector<std::int32_t> entities;
212 entities.reserve(int_ents_mesh.size() / 2);
213 if (codim == 0)
214 {
215 // In the codim 0 case, we need to map from cells in
216 // `this->mesh()` to cells in the argument/coefficient mesh, so
217 // here we extract the cells.
218 for (std::size_t i = 0; i < int_ents_mesh.size(); i += 2)
219 entities.push_back(int_ents_mesh[i]);
220 }
221 else if (codim == 1)
222 {
223 // In the codim 1 case, we need to map facets in `this->mesh()`
224 // to cells in the argument/coefficient mesh, so here we extract
225 // the facet index using the cell-to-facet connectivity.
226 for (std::size_t i = 0; i < int_ents_mesh.size(); i += 2)
227 {
228 entities.push_back(
229 c_to_f->links(int_ents_mesh[i])[int_ents_mesh[i + 1]]);
230 }
231 }
232 else
233 throw std::runtime_error("Codimension > 1 not supported.");
234
235 // Map from entity indices in `this->mesh()` to the corresponding
236 // cell indices in the argument/coefficient mesh
237 std::vector<std::int32_t> cells_mesh0
238 = emap.sub_topology_to_topology(entities, inverse);
239
240 // Create a list of (cell, local_facet_index) pairs in the
241 // argument/coefficient domain. Since `create_submesh`preserves
242 // the local facet index (with respect to the cell), we can use
243 // the local facet indices from the input integration entities
244 std::vector<std::int32_t> e = int_ents_mesh;
245 for (std::size_t i = 0; i < cells_mesh0.size(); ++i)
246 e[2 * i] = cells_mesh0[i];
247
248 return e;
249 };
250
251 for (auto& space : _function_spaces)
252 {
253 // Working map: [integral type, integral_idx, kernel_idx]->entities
254 std::map<std::tuple<IntegralType, int, int>,
255 std::variant<std::vector<std::int32_t>,
256 std::span<const std::int32_t>>>
257 vdata;
258
259 if (auto mesh0 = space->mesh(); mesh0 == _mesh)
260 {
261 for (auto& [key, integral] : _integrals)
262 vdata.insert({key, std::span(integral.entities)});
263 }
264 else
265 {
266 // Find correct entity map
267 const mesh::EntityMap& emap = get_entity_map(mesh0);
268
269 // Determine direction of the map. We need to map from
270 // `this->mesh()` to `mesh0`, so if `emap->sub_topology()` isn't
271 // the source topology, we need the inverse map
272 bool inverse = emap.sub_topology() == mesh0->topology();
273 for (auto& [key, itg] : _integrals)
274 {
275 auto [type, idx, kernel_idx] = key;
276 std::vector<std::int32_t> e;
277 if (type == IntegralType::cell)
278 e = emap.sub_topology_to_topology(itg.entities, inverse);
279 else if (type == IntegralType::exterior_facet
281 {
282 const mesh::Topology topology = *_mesh->topology();
283 int tdim = topology.dim();
284 assert(mesh0);
285 int codim = tdim - mesh0->topology()->dim();
286 assert(codim >= 0);
287 auto c_to_f = topology.connectivity(tdim, tdim - 1);
288 assert(c_to_f);
289
290 e = compute_facet_domains(itg.entities, codim, c_to_f, emap,
291 inverse);
292 }
293 else
294 throw std::runtime_error("Integral type not supported.");
295
296 vdata.insert({key, std::move(e)});
297 }
298 }
299
300 _edata.push_back(vdata);
301 }
302
303 for (auto& [key, integral] : _integrals)
304 {
305 auto [type, idx, kernel_idx] = key;
306 for (int c : integral.coeffs)
307 {
308 if (auto mesh0 = coefficients.at(c)->function_space()->mesh();
309 mesh0 == _mesh)
310 {
311 _cdata.insert({{type, idx, c}, std::span(integral.entities)});
312 }
313 else
314 {
315 // Find correct entity map and determine direction of the map
316 const mesh::EntityMap& emap = get_entity_map(mesh0);
317 bool inverse = emap.sub_topology() == mesh0->topology();
318
319 std::vector<std::int32_t> e;
320 if (type == IntegralType::cell)
321 e = emap.sub_topology_to_topology(integral.entities, inverse);
322 else if (type == IntegralType::exterior_facet
324 {
325 const mesh::Topology topology = *_mesh->topology();
326 int tdim = topology.dim();
327 assert(mesh0);
328 int codim = tdim - mesh0->topology()->dim();
329 auto c_to_f = topology.connectivity(tdim, tdim - 1);
330 assert(c_to_f);
331
332 e = compute_facet_domains(integral.entities, codim, c_to_f, emap,
333 inverse);
334 }
335 else
336 throw std::runtime_error("Integral type not supported.");
337 _cdata.insert({{type, idx, c}, std::move(e)});
338 }
339 }
340 }
341 }
342
344 Form(const Form& form) = delete;
345
347 Form(Form&& form) = default;
348
350 virtual ~Form() = default;
351
357 int rank() const { return _function_spaces.size(); }
358
361 std::shared_ptr<const mesh::Mesh<geometry_type>> mesh() const
362 {
363 return _mesh;
364 }
365
368 const std::vector<std::shared_ptr<const FunctionSpace<geometry_type>>>&
370 {
371 return _function_spaces;
372 }
373
383 std::function<void(scalar_type*, const scalar_type*, const scalar_type*,
384 const geometry_type*, const int*, const uint8_t*, void*)>
385 kernel(IntegralType type, int id, int kernel_idx) const
386 {
387 auto it = _integrals.find({type, id, kernel_idx});
388 if (it == _integrals.end())
389 throw std::runtime_error("Requested integral kernel not found.");
390 return it->second.kernel;
391 }
392
395 std::set<IntegralType> integral_types() const
396 {
397 std::vector<IntegralType> set_data;
398 std::ranges::transform(_integrals, std::back_inserter(set_data),
399 [](auto& x) { return std::get<0>(x.first); });
400 return std::set<IntegralType>(set_data.begin(), set_data.end());
401 }
402
413 std::vector<int> active_coeffs(IntegralType type, int id) const
414 {
415 auto it = std::ranges::find_if(_integrals,
416 [type, id](auto& x)
417 {
418 auto [t, id_, kernel_idx] = x.first;
419 return t == type and id_ == id;
420 });
421 if (it == _integrals.end())
422 throw std::runtime_error("Could not find active coefficient list.");
423 return it->second.coeffs;
424 }
425
439 int num_integrals(IntegralType type, int kernel_idx) const
440 {
441
442 int count = std::count_if(_integrals.begin(), _integrals.end(),
443 [type, kernel_idx](auto& x)
444 {
445 auto [t, id, k_idx] = x.first;
446 return t == type and k_idx == kernel_idx;
447 });
448 return count;
449 }
450
482 std::span<const std::int32_t> domain(IntegralType type, int idx,
483 int kernel_idx) const
484 {
485 auto it = _integrals.find({type, idx, kernel_idx});
486 if (it == _integrals.end())
487 throw std::runtime_error("Requested domain not found.");
488 return it->second.entities;
489 }
490
525 std::span<const std::int32_t> domain_arg(IntegralType type, int rank, int idx,
526 int kernel_idx) const
527 {
528 auto it = _edata.at(rank).find({type, idx, kernel_idx});
529 if (it == _edata.at(rank).end())
530 throw std::runtime_error("Requested domain for argument not found.");
531 try
532 {
533 return std::get<std::span<const std::int32_t>>(it->second);
534 }
535 catch (std::bad_variant_access& e)
536 {
537 return std::get<std::vector<std::int32_t>>(it->second);
538 }
539 }
540
555 std::span<const std::int32_t> domain_coeff(IntegralType type, int idx,
556 int c) const
557 {
558 auto it = _cdata.find({type, idx, c});
559 if (it == _cdata.end())
560 throw std::runtime_error("No domain for requested integral.");
561 try
562 {
563 return std::get<std::span<const std::int32_t>>(it->second);
564 }
565 catch (std::bad_variant_access& e)
566 {
567 return std::get<std::vector<std::int32_t>>(it->second);
568 }
569 }
570
572 const std::vector<
573 std::shared_ptr<const Function<scalar_type, geometry_type>>>&
575 {
576 return _coefficients;
577 }
578
582 bool needs_facet_permutations() const { return _needs_facet_permutations; }
583
588 std::vector<int> coefficient_offsets() const
589 {
590 std::vector<int> n{0};
591 for (auto& c : _coefficients)
592 {
593 if (!c)
594 throw std::runtime_error("Not all form coefficients have been set.");
595 n.push_back(n.back() + c->function_space()->element()->space_dimension());
596 }
597 return n;
598 }
599
601 const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
602 constants() const
603 {
604 return _constants;
605 }
606
607private:
608 // Function spaces (one for each argument)
609 std::vector<std::shared_ptr<const FunctionSpace<geometry_type>>>
610 _function_spaces;
611
612 // Integrals (integral type, id, celltype)
613 std::map<std::tuple<IntegralType, int, int>,
615 _integrals;
616
617 // The mesh
618 std::shared_ptr<const mesh::Mesh<geometry_type>> _mesh;
619
620 // Form coefficients
621 std::vector<std::shared_ptr<const Function<scalar_type, geometry_type>>>
622 _coefficients;
623
624 // Constants associated with the Form
625 std::vector<std::shared_ptr<const Constant<scalar_type>>> _constants;
626
627 // True if permutation data needs to be passed into these integrals
628 bool _needs_facet_permutations;
629
630 // Mapped domain index data for argument functions.
631 //
632 // Consider:
633 //
634 // entities = this->domain(IntegralType, integral(id), kernel_idx];
635 // entities0 = _edata[0][IntegralType, integral(id), coefficient_index];
636 //
637 // Then `entities[i]` is a mesh entity index (e.g., cell index) in
638 // `_mesh`, and `entities0[i]` is the index of the same entity but in
639 // the mesh associated with the argument 0 (test function) space.
640 std::vector<std::map<
641 std::tuple<IntegralType, int, int>,
642 std::variant<std::vector<std::int32_t>, std::span<const std::int32_t>>>>
643 _edata;
644
645 // Mapped domain index data for coefficient functions.
646 //
647 // Consider:
648 //
649 // entities = this->domain(IntegralType, integral(id), kernel_idx];
650 // entities0 = _cdata[IntegralType, integral(id), coefficient_index];
651 //
652 // Then `entities[i]` is a mesh entity index (e.g., cell index) in
653 // `_mesh`, and `entities0[i]` is the index of the same entity but in
654 // the mesh associated with the coefficient Function.
655 std::map<
656 std::tuple<IntegralType, int, int>,
657 std::variant<std::vector<std::int32_t>, std::span<const std::int32_t>>>
658 _cdata;
659};
660} // namespace dolfinx::fem
Constant (in space) value which can be attached to a Form.
Definition Constant.h:22
U geometry_type
Geometry type.
Definition Form.h:122
int num_integrals(IntegralType type, int kernel_idx) const
Get number of integrals (kernels) for a given integral type and kernel index.
Definition Form.h:439
const std::vector< std::shared_ptr< const Function< scalar_type, geometry_type > > > & coefficients() const
Access coefficients.
Definition Form.h:574
bool needs_facet_permutations() const
Get bool indicating whether permutation data needs to be passed into these integrals.
Definition Form.h:582
Form(Form &&form)=default
Move constructor.
int rank() const
Rank of the form.
Definition Form.h:357
Form(const Form &form)=delete
Copy constructor.
virtual ~Form()=default
Destructor.
Form(const std::vector< std::shared_ptr< const FunctionSpace< geometry_type > > > &V, X &&integrals, std::shared_ptr< const mesh::Mesh< geometry_type > > mesh, const std::vector< std::shared_ptr< const Function< scalar_type, geometry_type > > > &coefficients, const std::vector< std::shared_ptr< const Constant< scalar_type > > > &constants, bool needs_facet_permutations, const std::vector< std::reference_wrapper< const mesh::EntityMap > > &entity_maps)
Create a finite element form.
Definition Form.h:159
const std::vector< std::shared_ptr< const Constant< scalar_type > > > & constants() const
Access constants.
Definition Form.h:602
std::function< void(scalar_type *, const scalar_type *, const scalar_type *, const geometry_type *, const int *, const uint8_t *, void *)> kernel(IntegralType type, int id, int kernel_idx) const
Get the kernel function for an integral.
Definition Form.h:385
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell.
Definition Form.h:588
std::span< const std::int32_t > domain_coeff(IntegralType type, int idx, int c) const
Coefficient function mesh integration entity indices.
Definition Form.h:555
std::vector< int > active_coeffs(IntegralType type, int id) const
Indices of coefficients that are active for a given integral (kernel).
Definition Form.h:413
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
Common mesh for the form (the 'integration domain').
Definition Form.h:361
std::set< IntegralType > integral_types() const
Get types of integrals in the form.
Definition Form.h:395
const std::vector< std::shared_ptr< const FunctionSpace< geometry_type > > > & function_spaces() const
Function spaces for all arguments.
Definition Form.h:369
std::span< const std::int32_t > domain_arg(IntegralType type, int rank, int idx, int kernel_idx) const
Argument function mesh integration entity indices.
Definition Form.h:525
std::span< const std::int32_t > domain(IntegralType type, int idx, int kernel_idx) const
Mesh entity indices to integrate over for a given integral (kernel).
Definition Form.h:482
T scalar_type
Scalar type.
Definition Form.h:119
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
Definition Function.h:47
A bidirectional map relating entities in one topology to another.
Definition EntityMap.h:20
std::vector< std::int32_t > sub_topology_to_topology(std::span< const std::int32_t > entities, bool inverse) const
Map entities between the sub-topology and the parent topology.
Definition EntityMap.cpp:30
std::shared_ptr< const Topology > sub_topology() const
Get the sub-topology.
Definition EntityMap.cpp:24
std::shared_ptr< const Topology > topology() const
Get the (parent) topology.
Definition EntityMap.cpp:19
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:41
std::shared_ptr< const graph::AdjacencyList< std::int32_t > > connectivity(std::array< int, 2 > d0, std::array< int, 2 > d1) const
Get the connectivity from entities of topological dimension d0 to dimension d1.
Definition Topology.cpp:840
int dim() const noexcept
Topological dimension of the mesh.
Definition Topology.cpp:786
Finite element method functionality.
Definition assemble_expression_impl.h:23
@ inverse
Inverse.
Definition FiniteElement.h:29
IntegralType
Type of integral.
Definition Form.h:39
@ vertex
Vertex.
Definition Form.h:43
@ interior_facet
Interior facet.
Definition Form.h:42
@ cell
Cell.
Definition Form.h:40
@ exterior_facet
Exterior facet.
Definition Form.h:41
Represents integral data, containing the kernel, and a list of entities to integrate over and the ind...
Definition Form.h:51
std::function< void(T *, const T *, const T *, const U *, const int *, const uint8_t *, void *)> kernel
The integration kernel.
Definition Form.h:75
integral_data(K &&kernel, V &&entities, W &&coeffs)
Create a structure to hold integral data.
Definition Form.h:66
std::vector< int > coeffs
Indices of coefficients (from the form) that are in this integral.
Definition Form.h:83
std::vector< std::int32_t > entities
The entities to integrate over for this integral. These are the entities in 'full' mesh.
Definition Form.h:79