10#include "FunctionSpace.h"
13#include <basix/mdspan.hpp>
16#include <dolfinx/common/IndexMap.h>
17#include <dolfinx/common/types.h>
18#include <dolfinx/mesh/EntityMap.h>
19#include <dolfinx/mesh/Mesh.h>
32template <dolfinx::scalar T>
34template <dolfinx::scalar T, std::
floating_po
int U>
49template <dolfinx::scalar T, std::
floating_po
int U = scalar_value_t<T>>
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>,
73 std::function<void(T*,
const T*,
const T*,
const U*,
const int*,
74 const uint8_t*,
void*)>
114template <dolfinx::scalar T, std::
floating_po
int U = dolfinx::scalar_value_t<T>>
154 template <
typename X>
155 requires std::is_convertible_v<
156 std::remove_cvref_t<X>,
157 std::map<std::tuple<IntegralType, int, int>,
168 const std::vector<std::reference_wrapper<const mesh::EntityMap>>&
170 : _function_spaces(V), _integrals(std::forward<X>(integrals)),
175 throw std::runtime_error(
"Form Mesh is null.");
181 auto it = std::ranges::find_if(
185 return ((em.
topology() == mesh0->topology()
191 if (it == entity_maps.end())
193 throw std::runtime_error(
194 "Incompatible mesh. argument entity_maps must be provided.");
202 auto compute_facet_domains
203 = [&](
const auto& int_ents_mesh,
int codim,
const auto& c_to_f,
204 const auto& emap,
bool inverse)
211 std::vector<std::int32_t> entities;
212 entities.reserve(int_ents_mesh.size() / 2);
218 for (std::size_t i = 0; i < int_ents_mesh.size(); i += 2)
219 entities.push_back(int_ents_mesh[i]);
226 for (std::size_t i = 0; i < int_ents_mesh.size(); i += 2)
229 c_to_f->links(int_ents_mesh[i])[int_ents_mesh[i + 1]]);
233 throw std::runtime_error(
"Codimension > 1 not supported.");
237 std::vector<std::int32_t> cells_mesh0
238 = emap.sub_topology_to_topology(entities,
inverse);
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];
251 for (
auto& space : _function_spaces)
254 std::map<std::tuple<IntegralType, int, int>,
255 std::variant<std::vector<std::int32_t>,
256 std::span<const std::int32_t>>>
259 if (
auto mesh0 = space->mesh(); mesh0 == _mesh)
261 for (
auto& [key, integral] : _integrals)
262 vdata.insert({key, std::span(integral.entities)});
273 for (
auto& [key, itg] : _integrals)
275 auto [type, idx, kernel_idx] = key;
276 std::vector<std::int32_t> e;
283 int tdim = topology.
dim();
285 int codim = tdim - mesh0->topology()->dim();
290 e = compute_facet_domains(itg.entities, codim, c_to_f, emap,
294 throw std::runtime_error(
"Integral type not supported.");
296 vdata.insert({key, std::move(e)});
300 _edata.push_back(vdata);
303 for (
auto& [key, integral] : _integrals)
305 auto [type, idx, kernel_idx] = key;
306 for (
int c : integral.coeffs)
308 if (
auto mesh0 =
coefficients.at(c)->function_space()->mesh();
311 _cdata.insert({{type, idx, c}, std::span(integral.entities)});
319 std::vector<std::int32_t> e;
326 int tdim = topology.
dim();
328 int codim = tdim - mesh0->topology()->dim();
332 e = compute_facet_domains(integral.entities, codim, c_to_f, emap,
336 throw std::runtime_error(
"Integral type not supported.");
337 _cdata.insert({{type, idx, c}, std::move(e)});
357 int rank()
const {
return _function_spaces.size(); }
361 std::shared_ptr<const mesh::Mesh<geometry_type>>
mesh()
const
368 const std::vector<std::shared_ptr<const FunctionSpace<geometry_type>>>&
371 return _function_spaces;
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;
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());
415 auto it = std::ranges::find_if(_integrals,
418 auto [t, id_, kernel_idx] = x.first;
419 return t == type and id_ == id;
421 if (it == _integrals.end())
422 throw std::runtime_error(
"Could not find active coefficient list.");
423 return it->second.coeffs;
442 int count = std::count_if(_integrals.begin(), _integrals.end(),
443 [type, kernel_idx](
auto& x)
445 auto [t, id, k_idx] = x.first;
446 return t == type and k_idx == kernel_idx;
483 int kernel_idx)
const
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;
526 int kernel_idx)
const
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.");
533 return std::get<std::span<const std::int32_t>>(it->second);
535 catch (std::bad_variant_access& e)
537 return std::get<std::vector<std::int32_t>>(it->second);
558 auto it = _cdata.find({type, idx, c});
559 if (it == _cdata.end())
560 throw std::runtime_error(
"No domain for requested integral.");
563 return std::get<std::span<const std::int32_t>>(it->second);
565 catch (std::bad_variant_access& e)
567 return std::get<std::vector<std::int32_t>>(it->second);
573 std::shared_ptr<const Function<scalar_type, geometry_type>>>&
576 return _coefficients;
590 std::vector<int> n{0};
591 for (
auto& c : _coefficients)
594 throw std::runtime_error(
"Not all form coefficients have been set.");
595 n.push_back(n.back() + c->function_space()->element()->space_dimension());
601 const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
609 std::vector<std::shared_ptr<const FunctionSpace<geometry_type>>>
613 std::map<std::tuple<IntegralType, int, int>,
618 std::shared_ptr<const mesh::Mesh<geometry_type>> _mesh;
621 std::vector<std::shared_ptr<const Function<scalar_type, geometry_type>>>
625 std::vector<std::shared_ptr<const Constant<scalar_type>>> _constants;
628 bool _needs_facet_permutations;
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>>>>
656 std::tuple<IntegralType, int, int>,
657 std::variant<std::vector<std::int32_t>, std::span<const std::int32_t>>>
Constant (in space) value which can be attached to a Form.
Definition Constant.h:22
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
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