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>
214template <dolfinx::scalar T, std::
floating_po
int U>
216 std::span<T> b,
const Form<T, U>& L, std::span<const T> constants,
217 const std::map<std::pair<IntegralType, int>,
218 std::pair<std::span<const T>,
int>>& coefficients)
220 impl::assemble_vector(b, L, constants, coefficients);
227template <dolfinx::scalar T, std::
floating_po
int U>
314template <dolfinx::scalar T, std::
floating_po
int U>
317 std::vector<std::optional<std::reference_wrapper<
const Form<T, U>>>> a,
318 const std::vector<std::span<const T>>& constants,
319 const std::vector<std::map<std::pair<IntegralType, int>,
320 std::pair<std::span<const T>,
int>>>& coeffs,
323 const std::vector<std::span<const T>>& x0, T alpha)
326 if (std::ranges::all_of(a, [](
auto ai) {
return !ai; }))
329 impl::apply_lifting<T>(b, a, constants, coeffs, bcs1, x0, alpha);
360template <dolfinx::scalar T, std::
floating_po
int U>
363 std::vector<std::optional<std::reference_wrapper<
const Form<T, U>>>> a,
366 const std::vector<std::span<const T>>& x0, T alpha)
369 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>,
int>>>
371 std::vector<std::vector<T>> constants;
378 coeffs.push_back(coefficients);
383 coeffs.emplace_back();
384 constants.emplace_back();
388 std::vector<std::span<const T>> _constants(constants.begin(),
390 std::vector<std::map<std::pair<IntegralType, int>,
391 std::pair<std::span<const T>,
int>>>
393 std::ranges::transform(coeffs, std::back_inserter(_coeffs),
413template <dolfinx::scalar T, std::
floating_po
int U>
416 std::span<const T> constants,
417 const std::map<std::pair<IntegralType, int>,
418 std::pair<std::span<const T>,
int>>& coefficients,
419 std::span<const std::int8_t> dof_marker0,
420 std::span<const std::int8_t> dof_marker1)
424 = md::mdspan<const scalar_value_t<T>,
425 md::extents<std::size_t, md::dynamic_extent, 3>>;
427 std::shared_ptr<const mesh::Mesh<U>>
mesh = a.mesh();
429 std::span x =
mesh->geometry().x();
430 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
432 impl::assemble_matrix(mat_add, a, mdspanx3_t(x.data(), x.size() / 3, 3),
433 constants, coefficients, dof_marker0, dof_marker1);
437 std::vector<scalar_value_t<T>> _x(x.begin(), x.end());
438 impl::assemble_matrix(mat_add, a, mdspanx3_t(_x.data(), _x.size() / 3, 3),
439 constants, coefficients, dof_marker0, dof_marker1);
450template <dolfinx::scalar T, std::
floating_po
int U>
452 auto mat_add,
const Form<T, U>& a, std::span<const T> constants,
453 const std::map<std::pair<IntegralType, int>,
454 std::pair<std::span<const T>,
int>>& coefficients,
460 auto map0 = a.function_spaces().at(0)->dofmaps(0)->index_map;
461 auto map1 = a.function_spaces().at(1)->dofmaps(0)->index_map;
462 auto bs0 = a.function_spaces().at(0)->dofmaps(0)->index_map_bs();
463 auto bs1 = a.function_spaces().at(1)->dofmaps(0)->index_map_bs();
466 std::vector<std::int8_t> dof_marker0, dof_marker1;
468 std::int32_t dim0 = bs0 * (map0->size_local() + map0->num_ghosts());
470 std::int32_t dim1 = bs1 * (map1->size_local() + map1->num_ghosts());
471 for (std::size_t k = 0; k < bcs.size(); ++k)
473 assert(bcs[k].get().function_space());
474 if (a.function_spaces().at(0)->contains(*bcs[k].get().function_space()))
476 dof_marker0.resize(dim0,
false);
477 bcs[k].get().mark_dofs(dof_marker0);
480 if (a.function_spaces().at(1)->contains(*bcs[k].get().function_space()))
482 dof_marker1.resize(dim1,
false);
483 bcs[k].get().mark_dofs(dof_marker1);
497template <dolfinx::scalar T, std::
floating_po
int U>
523template <dolfinx::scalar T, std::
floating_po
int U>
525 std::span<const std::int8_t> dof_marker0,
526 std::span<const std::int8_t> dof_marker1)
552template <dolfinx::scalar T>
556 for (std::size_t i = 0; i < rows.size(); ++i)
558 std::span diag_span(&diagonal, 1);
559 set_fn(rows.subspan(i, 1), rows.subspan(i, 1), diag_span);
579template <dolfinx::scalar T, std::
floating_po
int U>
587 if (V.contains(*bc.get().function_space()))
589 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:414
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:553
void assemble_vector(std::span< T > 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:215
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
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
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
void apply_lifting(std::span< T > 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:315
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
Functions supporting the packing of coefficient data.