10#include "FunctionSpace.h"
11#include "assemble_expression_impl.h"
12#include "assemble_matrix_impl.h"
13#include "assemble_scalar_impl.h"
14#include "assemble_vector_impl.h"
19#include <basix/mdspan.hpp>
21#include <dolfinx/common/types.h>
33template <dolfinx::scalar T, std::
floating_po
int U>
35template <dolfinx::scalar T, std::
floating_po
int U>
37template <dolfinx::scalar T, std::
floating_po
int U>
39template <std::
floating_po
int T>
63template <dolfinx::scalar T, std::
floating_po
int U>
66 md::mdspan<
const T, md::dextents<std::size_t, 2>> coeffs,
73 auto [X, Xshape] = e.X();
74 impl::tabulate_expression(values, e.kernel(), Xshape, e.value_size(), coeffs,
75 constants,
mesh, entities, element);
94template <dolfinx::scalar T, std::
floating_po
int U>
99 std::pair<std::reference_wrapper<const FiniteElement<U>>, std::size_t>>
100 element = std::nullopt;
101 if (
auto V = e.argument_space(); V)
103 std::size_t num_argument_dofs
104 = V->dofmap()->element_dof_layout().num_dofs() * V->dofmap()->bs();
105 assert(V->element());
106 element = {std::cref(*V->element()), num_argument_dofs};
109 std::vector<int> coffsets = e.coefficient_offsets();
110 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients
112 std::vector<T> coeffs(entities.extent(0) * coffsets.back());
113 int cstride = coffsets.back();
115 std::vector<std::reference_wrapper<const Function<T, U>>> c;
116 std::ranges::transform(coefficients, std::back_inserter(c),
123 values, e, md::mdspan(coeffs.data(), entities.size(), cstride),
124 std::span<const T>(constants),
mesh, entities, element);
130template <dolfinx::scalar T>
131std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>,
int>>
133 std::pair<std::vector<T>,
int>>& coeffs)
135 using Key =
typename std::remove_reference_t<
decltype(coeffs)>::key_type;
136 std::map<Key, std::pair<std::span<const T>,
int>> c;
137 std::ranges::transform(
138 coeffs, std::inserter(c, c.end()),
139 [](
auto& e) ->
typename decltype(c)::value_type
140 { return {e.first, {e.second.first, e.second.second}}; });
157template <dolfinx::scalar T, std::
floating_po
int U>
159 const Form<T, U>& M, std::span<const T> constants,
160 const std::map<std::pair<IntegralType, int>,
161 std::pair<std::span<const T>,
int>>& coefficients)
164 = md::mdspan<const scalar_value_t<T>,
165 md::extents<std::size_t, md::dynamic_extent, 3>>;
167 std::shared_ptr<const mesh::Mesh<U>>
mesh = M.mesh();
169 std::span x =
mesh->geometry().x();
170 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
172 return impl::assemble_scalar(M,
mesh->geometry().dofmap(),
173 mdspanx3_t(x.data(), x.size() / 3, 3),
174 constants, coefficients);
178 std::vector<scalar_value_t<T>> _x(x.begin(), x.end());
179 return impl::assemble_scalar(M,
mesh->geometry().dofmap(),
180 mdspanx3_t(_x.data(), _x.size() / 3, 3),
181 constants, coefficients);
192template <dolfinx::scalar T, std::
floating_po
int U>
215template <
typename V, std::floating_point U,
217 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
219 V&& b,
const Form<T, U>& L, std::span<const T> constants,
220 const std::map<std::pair<IntegralType, int>,
221 std::pair<std::span<const T>,
int>>& coefficients)
223 impl::assemble_vector(b, L, constants, coefficients);
232template <
typename V, std::floating_point U,
234 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
322 std::floating_point U
323 = scalar_value_t<typename std::remove_cvref_t<V>::value_type>,
325 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
328 std::vector<std::optional<std::reference_wrapper<
const Form<T, U>>>> a,
329 const std::vector<std::span<const T>>& constants,
330 const std::vector<std::map<std::pair<IntegralType, int>,
331 std::pair<std::span<const T>,
int>>>& coeffs,
334 const std::vector<std::span<const T>>& x0, T alpha)
337 if (std::ranges::all_of(a, [](
auto ai) {
return !ai; }))
340 impl::apply_lifting(b, a, constants, coeffs, bcs1, x0, alpha);
372 std::floating_point U
373 = scalar_value_t<typename std::remove_cvref_t<V>::value_type>,
375 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
378 std::vector<std::optional<std::reference_wrapper<
const Form<T, U>>>> a,
381 const std::vector<std::span<const T>>& x0, T alpha)
384 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>,
int>>>
386 std::vector<std::vector<T>> constants;
393 coeffs.push_back(coefficients);
398 coeffs.emplace_back();
399 constants.emplace_back();
403 std::vector<std::span<const T>> _constants(constants.begin(),
405 std::vector<std::map<std::pair<IntegralType, int>,
406 std::pair<std::span<const T>,
int>>>
408 std::ranges::transform(coeffs, std::back_inserter(_coeffs),
428template <dolfinx::scalar T, std::
floating_po
int U>
431 std::span<const T> constants,
432 const std::map<std::pair<IntegralType, int>,
433 std::pair<std::span<const T>,
int>>& coefficients,
434 std::span<const std::int8_t> dof_marker0,
435 std::span<const std::int8_t> dof_marker1)
439 = md::mdspan<const scalar_value_t<T>,
440 md::extents<std::size_t, md::dynamic_extent, 3>>;
442 std::shared_ptr<const mesh::Mesh<U>>
mesh = a.mesh();
444 std::span x =
mesh->geometry().x();
445 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
447 impl::assemble_matrix(mat_add, a, mdspanx3_t(x.data(), x.size() / 3, 3),
448 constants, coefficients, dof_marker0, dof_marker1);
452 std::vector<scalar_value_t<T>> _x(x.begin(), x.end());
453 impl::assemble_matrix(mat_add, a, mdspanx3_t(_x.data(), _x.size() / 3, 3),
454 constants, coefficients, dof_marker0, dof_marker1);
465template <dolfinx::scalar T, std::
floating_po
int U>
467 auto mat_add,
const Form<T, U>& a, std::span<const T> constants,
468 const std::map<std::pair<IntegralType, int>,
469 std::pair<std::span<const T>,
int>>& coefficients,
475 auto map0 = a.function_spaces().at(0)->dofmaps(0)->index_map;
476 auto map1 = a.function_spaces().at(1)->dofmaps(0)->index_map;
477 auto bs0 = a.function_spaces().at(0)->dofmaps(0)->index_map_bs();
478 auto bs1 = a.function_spaces().at(1)->dofmaps(0)->index_map_bs();
481 std::vector<std::int8_t> dof_marker0, dof_marker1;
483 std::int32_t dim0 = bs0 * (map0->size_local() + map0->num_ghosts());
485 std::int32_t dim1 = bs1 * (map1->size_local() + map1->num_ghosts());
486 for (std::size_t k = 0; k < bcs.size(); ++k)
488 assert(bcs[k].get().function_space());
489 if (a.function_spaces().at(0)->contains(*bcs[k].get().function_space()))
491 dof_marker0.resize(dim0,
false);
492 bcs[k].get().mark_dofs(dof_marker0);
495 if (a.function_spaces().at(1)->contains(*bcs[k].get().function_space()))
497 dof_marker1.resize(dim1,
false);
498 bcs[k].get().mark_dofs(dof_marker1);
512template <dolfinx::scalar T, std::
floating_po
int U>
538template <dolfinx::scalar T, std::
floating_po
int U>
540 std::span<const std::int8_t> dof_marker0,
541 std::span<const std::int8_t> dof_marker1)
567template <dolfinx::scalar T>
571 for (std::size_t i = 0; i < rows.size(); ++i)
573 std::span diag_span(&diagonal, 1);
574 set_fn(rows.subspan(i, 1), rows.subspan(i, 1), diag_span);
594template <dolfinx::scalar T, std::
floating_po
int U>
602 if (V.contains(*bc.get().function_space()))
604 const auto [dofs, range] = bc.get().dof_indices();
Definition DirichletBC.h:255
An Expression represents a mathematical expression evaluated at a pre-defined points on a reference c...
Definition Expression.h:41
Model of a finite element.
Definition FiniteElement.h:57
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Concept for mdspan of rank 1 or 2.
Definition traits.h:36
Matrix accumulate/set concept for functions that can be used in assemblers to accumulate or set value...
Definition utils.h:28
Functions supporting finite element method operations.
Finite element method functionality.
Definition assemble_expression_impl.h:23
std::pair< std::vector< T >, int > allocate_coefficient_storage(const Form< T, U > &form, IntegralType integral_type, int id)
Allocate storage for coefficients of a pair (integral_type, / id) from a Form.
Definition pack.h:170
void assemble_matrix(la::MatSet< T > auto mat_add, const Form< T, U > &a, std::span< const T > constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients, std::span< const std::int8_t > dof_marker0, std::span< const std::int8_t > dof_marker1)
Assemble bilinear form into a matrix. Matrix must already be initialised. Does not zero or finalise t...
Definition assembler.h:429
T assemble_scalar(const Form< T, U > &M, std::span< const T > constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients)
Assemble functional into scalar.
Definition assembler.h:158
void set_diagonal(auto set_fn, std::span< const std::int32_t > rows, T diagonal=1.0)
Sets a value to the diagonal of a matrix for specified rows.
Definition assembler.h:568
void tabulate_expression(std::span< T > values, const fem::Expression< T, U > &e, md::mdspan< const T, md::dextents< std::size_t, 2 > > coeffs, std::span< const T > constants, const mesh::Mesh< U > &mesh, fem::MDSpan2 auto entities, std::optional< std::pair< std::reference_wrapper< const FiniteElement< U > >, std::size_t > > element)
Evaluate an Expression on cells or facets.
Definition assembler.h:64
void apply_lifting(V &&b, std::vector< std::optional< std::reference_wrapper< const Form< T, U > > > > a, const std::vector< std::span< const T > > &constants, const std::vector< std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > > &coeffs, const std::vector< std::vector< std::reference_wrapper< const DirichletBC< T, U > > > > &bcs1, const std::vector< std::span< const T > > &x0, T alpha)
Modify the right-hand side vector to account for constraints (Dirichlet boundary condition constraint...
Definition assembler.h:326
std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > make_coefficients_span(const std::map< std::pair< IntegralType, int >, std::pair< std::vector< T >, int > > &coeffs)
Create a map of std::spans from a map of std::vectors.
Definition assembler.h:132
void pack_coefficients(const Form< T, U > &form, std::map< std::pair< IntegralType, int >, std::pair< std::vector< T >, int > > &coeffs)
Pack coefficients of a Form.
Definition pack.h:225
void assemble_vector(V &&b, const Form< T, U > &L, std::span< const T > constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients)
Assemble linear form into a vector.
Definition assembler.h:218
std::vector< T > pack_constants(std::vector< std::reference_wrapper< const fem::Constant< T > > > c)
Pack constants of an Expression or Form into a single array ready for assembly.
Definition pack.h:379
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
Functions supporting the packing of coefficient data.