DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
assembler.h
Go to the documentation of this file.
1// Copyright (C) 2018-2025 Garth N. Wells
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 "Function.h"
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"
15#include "pack.h"
16#include "traits.h"
17#include "utils.h"
18#include <algorithm>
19#include <basix/mdspan.hpp>
20#include <cstdint>
21#include <dolfinx/common/types.h>
22#include <memory>
23#include <optional>
24#include <span>
25#include <vector>
26
30
31namespace dolfinx::fem
32{
33template <dolfinx::scalar T, std::floating_point U>
34class DirichletBC;
35template <dolfinx::scalar T, std::floating_point U>
36class Expression;
37template <dolfinx::scalar T, std::floating_point U>
38class Form;
39template <std::floating_point T>
40class FunctionSpace;
41
63template <dolfinx::scalar T, std::floating_point U>
65 std::span<T> values, const fem::Expression<T, U>& e,
66 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
67 std::span<const T> constants, const mesh::Mesh<U>& mesh,
68 fem::MDSpan2 auto entities,
69 std::optional<
70 std::pair<std::reference_wrapper<const FiniteElement<U>>, std::size_t>>
71 element)
72{
73 auto [X, Xshape] = e.X();
74 impl::tabulate_expression(values, e.kernel(), Xshape, e.value_size(), coeffs,
75 constants, mesh, entities, element);
76}
77
94template <dolfinx::scalar T, std::floating_point U>
95void tabulate_expression(std::span<T> values, const fem::Expression<T, U>& e,
96 const mesh::Mesh<U>& mesh, fem::MDSpan2 auto entities)
97{
98 std::optional<
99 std::pair<std::reference_wrapper<const FiniteElement<U>>, std::size_t>>
100 element = std::nullopt;
101 if (auto V = e.argument_space(); V)
102 {
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};
107 }
108
109 std::vector<int> coffsets = e.coefficient_offsets();
110 const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients
111 = e.coefficients();
112 std::vector<T> coeffs(entities.extent(0) * coffsets.back());
113 int cstride = coffsets.back();
114 {
115 std::vector<std::reference_wrapper<const Function<T, U>>> c;
116 std::ranges::transform(coefficients, std::back_inserter(c),
117 [](auto c) -> const Function<T, U>& { return *c; });
118 fem::pack_coefficients(c, coffsets, entities, std::span(coeffs));
119 }
120 std::vector<T> constants = fem::pack_constants(e);
121
123 values, e, md::mdspan(coeffs.data(), entities.size(), cstride),
124 std::span<const T>(constants), mesh, entities, element);
125}
126
127// -- Helper functions -----------------------------------------------------
128
130template <dolfinx::scalar T>
131std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>>
132make_coefficients_span(const std::map<std::pair<IntegralType, int>,
133 std::pair<std::vector<T>, int>>& coeffs)
134{
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}}; });
141 return c;
142}
143
144// -- Scalar ----------------------------------------------------------------
145
157template <dolfinx::scalar T, std::floating_point 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)
162{
163 using mdspanx3_t
164 = md::mdspan<const scalar_value_t<T>,
165 md::extents<std::size_t, md::dynamic_extent, 3>>;
166
167 std::shared_ptr<const mesh::Mesh<U>> mesh = M.mesh();
168 assert(mesh);
169 std::span x = mesh->geometry().x();
170 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
171 {
172 return impl::assemble_scalar(M, mesh->geometry().dofmap(),
173 mdspanx3_t(x.data(), x.size() / 3, 3),
174 constants, coefficients);
175 }
176 else
177 {
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);
182 }
183}
184
192template <dolfinx::scalar T, std::floating_point U>
194{
195 const std::vector<T> constants = pack_constants(M);
196 auto coefficients = allocate_coefficient_storage(M);
197 pack_coefficients(M, coefficients);
198 return assemble_scalar(M, std::span(constants),
199 make_coefficients_span(coefficients));
200}
201
202// -- Vectors ----------------------------------------------------------------
203
214// template <dolfinx::scalar T, std::floating_point U>
215template <typename V, std::floating_point U,
216 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
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)
222{
223 impl::assemble_vector(b, L, constants, coefficients);
224}
225
230// template <dolfinx::scalar T, std::floating_point U>
231// void assemble_vector(std::span<T> b, const Form<T, U>& L)
232template <typename V, std::floating_point U,
233 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
234 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
235void assemble_vector(V&& b, const Form<T, U>& L)
236{
237 auto coefficients = allocate_coefficient_storage(L);
238 pack_coefficients(L, coefficients);
239 const std::vector<T> constants = pack_constants(L);
240 assemble_vector(b, L, std::span(constants),
241 make_coefficients_span(coefficients));
242}
243
321template <typename V,
322 std::floating_point U
323 = scalar_value_t<typename std::remove_cvref_t<V>::value_type>,
324 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
325 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
327 V&& b,
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,
332 const std::vector<
333 std::vector<std::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
334 const std::vector<std::span<const T>>& x0, T alpha)
335{
336 // If all forms are null, there is nothing to do
337 if (std::ranges::all_of(a, [](auto ai) { return !ai; }))
338 return;
339
340 impl::apply_lifting(b, a, constants, coeffs, bcs1, x0, alpha);
341}
342
371template <typename V,
372 std::floating_point U
373 = scalar_value_t<typename std::remove_cvref_t<V>::value_type>,
374 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
375 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
377 V&& b,
378 std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> a,
379 const std::vector<
380 std::vector<std::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
381 const std::vector<std::span<const T>>& x0, T alpha)
382{
383 std::vector<
384 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>>
385 coeffs;
386 std::vector<std::vector<T>> constants;
387 for (auto _a : a)
388 {
389 if (_a)
390 {
391 auto coefficients = allocate_coefficient_storage(_a->get());
392 pack_coefficients(_a->get(), coefficients);
393 coeffs.push_back(coefficients);
394 constants.push_back(pack_constants(_a->get()));
395 }
396 else
397 {
398 coeffs.emplace_back();
399 constants.emplace_back();
400 }
401 }
402
403 std::vector<std::span<const T>> _constants(constants.begin(),
404 constants.end());
405 std::vector<std::map<std::pair<IntegralType, int>,
406 std::pair<std::span<const T>, int>>>
407 _coeffs;
408 std::ranges::transform(coeffs, std::back_inserter(_coeffs),
409 [](auto& c) { return make_coefficients_span(c); });
410
411 apply_lifting(b, a, _constants, _coeffs, bcs1, x0, alpha);
412}
413
414// -- Matrices ---------------------------------------------------------------
415
428template <dolfinx::scalar T, std::floating_point U>
430 la::MatSet<T> auto mat_add, const Form<T, U>& a,
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)
436
437{
438 using mdspanx3_t
439 = md::mdspan<const scalar_value_t<T>,
440 md::extents<std::size_t, md::dynamic_extent, 3>>;
441
442 std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
443 assert(mesh);
444 std::span x = mesh->geometry().x();
445 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
446 {
447 impl::assemble_matrix(mat_add, a, mdspanx3_t(x.data(), x.size() / 3, 3),
448 constants, coefficients, dof_marker0, dof_marker1);
449 }
450 else
451 {
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);
455 }
456}
457
465template <dolfinx::scalar T, std::floating_point 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,
470 const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs)
471{
472 // Index maps for dof ranges
473 // NOTE: For mixed-topology meshes, there will be multiple DOF maps,
474 // but the index maps are the same.
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();
479
480 // Build dof markers
481 std::vector<std::int8_t> dof_marker0, dof_marker1;
482 assert(map0);
483 std::int32_t dim0 = bs0 * (map0->size_local() + map0->num_ghosts());
484 assert(map1);
485 std::int32_t dim1 = bs1 * (map1->size_local() + map1->num_ghosts());
486 for (std::size_t k = 0; k < bcs.size(); ++k)
487 {
488 assert(bcs[k].get().function_space());
489 if (a.function_spaces().at(0)->contains(*bcs[k].get().function_space()))
490 {
491 dof_marker0.resize(dim0, false);
492 bcs[k].get().mark_dofs(dof_marker0);
493 }
494
495 if (a.function_spaces().at(1)->contains(*bcs[k].get().function_space()))
496 {
497 dof_marker1.resize(dim1, false);
498 bcs[k].get().mark_dofs(dof_marker1);
499 }
500 }
501
502 // Assemble
503 assemble_matrix(mat_add, a, constants, coefficients, dof_marker0,
504 dof_marker1);
505}
506
512template <dolfinx::scalar T, std::floating_point U>
514 auto mat_add, const Form<T, U>& a,
515 const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs)
516{
517 // Prepare constants and coefficients
518 const std::vector<T> constants = pack_constants(a);
519 auto coefficients = allocate_coefficient_storage(a);
520 pack_coefficients(a, coefficients);
521
522 // Assemble
523 assemble_matrix(mat_add, a, std::span(constants),
524 make_coefficients_span(coefficients), bcs);
525}
526
538template <dolfinx::scalar T, std::floating_point U>
539void assemble_matrix(auto mat_add, const Form<T, U>& a,
540 std::span<const std::int8_t> dof_marker0,
541 std::span<const std::int8_t> dof_marker1)
542
543{
544 // Prepare constants and coefficients
545 const std::vector<T> constants = pack_constants(a);
546 auto coefficients = allocate_coefficient_storage(a);
547 pack_coefficients(a, coefficients);
548
549 // Assemble
550 assemble_matrix(mat_add, a, std::span(constants),
551 make_coefficients_span(coefficients), dof_marker0,
552 dof_marker1);
553}
554
567template <dolfinx::scalar T>
568void set_diagonal(auto set_fn, std::span<const std::int32_t> rows,
569 T diagonal = 1.0)
570{
571 for (std::size_t i = 0; i < rows.size(); ++i)
572 {
573 std::span diag_span(&diagonal, 1);
574 set_fn(rows.subspan(i, 1), rows.subspan(i, 1), diag_span);
575 }
576}
577
594template <dolfinx::scalar T, std::floating_point U>
596 auto set_fn, const FunctionSpace<U>& V,
597 const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs,
598 T diagonal = 1.0)
599{
600 for (auto& bc : bcs)
601 {
602 if (V.contains(*bc.get().function_space()))
603 {
604 const auto [dofs, range] = bc.get().dof_indices();
605 set_diagonal(set_fn, dofs.first(range), diagonal);
606 }
607 }
608}
609
610} // namespace dolfinx::fem
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
A representation of finite element variational forms.
Definition Form.h:116
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
Definition Function.h:47
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
Definition types.h:20
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.