Loading [MathJax]/extensions/tex2jax.js
DOLFINx 0.10.0.0
DOLFINx C++ interface
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Modules Pages Concepts
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
214template <dolfinx::scalar T, std::floating_point 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)
219{
220 impl::assemble_vector(b, L, constants, coefficients);
221}
222
227template <dolfinx::scalar T, std::floating_point U>
228void assemble_vector(std::span<T> b, const Form<T, U>& L)
229{
230 auto coefficients = allocate_coefficient_storage(L);
231 pack_coefficients(L, coefficients);
232 const std::vector<T> constants = pack_constants(L);
233 assemble_vector(b, L, std::span(constants),
234 make_coefficients_span(coefficients));
235}
236
314template <dolfinx::scalar T, std::floating_point U>
316 std::span<T> b,
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,
321 const std::vector<
322 std::vector<std::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
323 const std::vector<std::span<const T>>& x0, T alpha)
324{
325 // If all forms are null, there is nothing to do
326 if (std::ranges::all_of(a, [](auto ai) { return !ai; }))
327 return;
328
329 impl::apply_lifting<T>(b, a, constants, coeffs, bcs1, x0, alpha);
330}
331
360template <dolfinx::scalar T, std::floating_point U>
362 std::span<T> b,
363 std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> a,
364 const std::vector<
365 std::vector<std::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
366 const std::vector<std::span<const T>>& x0, T alpha)
367{
368 std::vector<
369 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>>
370 coeffs;
371 std::vector<std::vector<T>> constants;
372 for (auto _a : a)
373 {
374 if (_a)
375 {
376 auto coefficients = allocate_coefficient_storage(_a->get());
377 pack_coefficients(_a->get(), coefficients);
378 coeffs.push_back(coefficients);
379 constants.push_back(pack_constants(_a->get()));
380 }
381 else
382 {
383 coeffs.emplace_back();
384 constants.emplace_back();
385 }
386 }
387
388 std::vector<std::span<const T>> _constants(constants.begin(),
389 constants.end());
390 std::vector<std::map<std::pair<IntegralType, int>,
391 std::pair<std::span<const T>, int>>>
392 _coeffs;
393 std::ranges::transform(coeffs, std::back_inserter(_coeffs),
394 [](auto& c) { return make_coefficients_span(c); });
395
396 apply_lifting(b, a, _constants, _coeffs, bcs1, x0, alpha);
397}
398
399// -- Matrices ---------------------------------------------------------------
400
413template <dolfinx::scalar T, std::floating_point U>
415 la::MatSet<T> auto mat_add, const Form<T, U>& a,
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)
421
422{
423 using mdspanx3_t
424 = md::mdspan<const scalar_value_t<T>,
425 md::extents<std::size_t, md::dynamic_extent, 3>>;
426
427 std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
428 assert(mesh);
429 std::span x = mesh->geometry().x();
430 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
431 {
432 impl::assemble_matrix(mat_add, a, mdspanx3_t(x.data(), x.size() / 3, 3),
433 constants, coefficients, dof_marker0, dof_marker1);
434 }
435 else
436 {
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);
440 }
441}
442
450template <dolfinx::scalar T, std::floating_point 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,
455 const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs)
456{
457 // Index maps for dof ranges
458 // NOTE: For mixed-topology meshes, there will be multiple DOF maps,
459 // but the index maps are the same.
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();
464
465 // Build dof markers
466 std::vector<std::int8_t> dof_marker0, dof_marker1;
467 assert(map0);
468 std::int32_t dim0 = bs0 * (map0->size_local() + map0->num_ghosts());
469 assert(map1);
470 std::int32_t dim1 = bs1 * (map1->size_local() + map1->num_ghosts());
471 for (std::size_t k = 0; k < bcs.size(); ++k)
472 {
473 assert(bcs[k].get().function_space());
474 if (a.function_spaces().at(0)->contains(*bcs[k].get().function_space()))
475 {
476 dof_marker0.resize(dim0, false);
477 bcs[k].get().mark_dofs(dof_marker0);
478 }
479
480 if (a.function_spaces().at(1)->contains(*bcs[k].get().function_space()))
481 {
482 dof_marker1.resize(dim1, false);
483 bcs[k].get().mark_dofs(dof_marker1);
484 }
485 }
486
487 // Assemble
488 assemble_matrix(mat_add, a, constants, coefficients, dof_marker0,
489 dof_marker1);
490}
491
497template <dolfinx::scalar T, std::floating_point U>
499 auto mat_add, const Form<T, U>& a,
500 const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs)
501{
502 // Prepare constants and coefficients
503 const std::vector<T> constants = pack_constants(a);
504 auto coefficients = allocate_coefficient_storage(a);
505 pack_coefficients(a, coefficients);
506
507 // Assemble
508 assemble_matrix(mat_add, a, std::span(constants),
509 make_coefficients_span(coefficients), bcs);
510}
511
523template <dolfinx::scalar T, std::floating_point U>
524void assemble_matrix(auto mat_add, const Form<T, U>& a,
525 std::span<const std::int8_t> dof_marker0,
526 std::span<const std::int8_t> dof_marker1)
527
528{
529 // Prepare constants and coefficients
530 const std::vector<T> constants = pack_constants(a);
531 auto coefficients = allocate_coefficient_storage(a);
532 pack_coefficients(a, coefficients);
533
534 // Assemble
535 assemble_matrix(mat_add, a, std::span(constants),
536 make_coefficients_span(coefficients), dof_marker0,
537 dof_marker1);
538}
539
552template <dolfinx::scalar T>
553void set_diagonal(auto set_fn, std::span<const std::int32_t> rows,
554 T diagonal = 1.0)
555{
556 for (std::size_t i = 0; i < rows.size(); ++i)
557 {
558 std::span diag_span(&diagonal, 1);
559 set_fn(rows.subspan(i, 1), rows.subspan(i, 1), diag_span);
560 }
561}
562
579template <dolfinx::scalar T, std::floating_point U>
581 auto set_fn, const FunctionSpace<U>& V,
582 const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs,
583 T diagonal = 1.0)
584{
585 for (auto& bc : bcs)
586 {
587 if (V.contains(*bc.get().function_space()))
588 {
589 const auto [dofs, range] = bc.get().dof_indices();
590 set_diagonal(set_fn, dofs.first(range), diagonal);
591 }
592 }
593}
594
595} // 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:175
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
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.