DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
assembler.h
1// Copyright (C) 2018-2022 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 "assemble_matrix_impl.h"
10#include "assemble_scalar_impl.h"
11#include "assemble_vector_impl.h"
12#include "traits.h"
13#include "utils.h"
14#include <algorithm>
15#include <cstdint>
16#include <dolfinx/common/types.h>
17#include <memory>
18#include <optional>
19#include <span>
20#include <vector>
21
22namespace dolfinx::fem
23{
24template <dolfinx::scalar T, std::floating_point U>
25class DirichletBC;
26template <dolfinx::scalar T, std::floating_point U>
27class Form;
28template <std::floating_point T>
29class FunctionSpace;
30
31// -- Helper functions -----------------------------------------------------
32
34template <dolfinx::scalar T>
35std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>>
36make_coefficients_span(const std::map<std::pair<IntegralType, int>,
37 std::pair<std::vector<T>, int>>& coeffs)
38{
39 using Key = typename std::remove_reference_t<decltype(coeffs)>::key_type;
40 std::map<Key, std::pair<std::span<const T>, int>> c;
41 std::ranges::transform(
42 coeffs, std::inserter(c, c.end()),
43 [](auto& e) -> typename decltype(c)::value_type
44 { return {e.first, {e.second.first, e.second.second}}; });
45 return c;
46}
47
48// -- Scalar ----------------------------------------------------------------
49
61template <dolfinx::scalar T, std::floating_point U>
62T assemble_scalar(
63 const Form<T, U>& M, std::span<const T> constants,
64 const std::map<std::pair<IntegralType, int>,
65 std::pair<std::span<const T>, int>>& coefficients)
66{
67 std::shared_ptr<const mesh::Mesh<U>> mesh = M.mesh();
68 assert(mesh);
69 if constexpr (std::is_same_v<U, scalar_value_type_t<T>>)
70 {
71 return impl::assemble_scalar(M, mesh->geometry().dofmap(),
72 mesh->geometry().x(), constants, coefficients);
73 }
74 else
75 {
76 auto x = mesh->geometry().x();
77 std::vector<scalar_value_type_t<T>> _x(x.begin(), x.end());
78 return impl::assemble_scalar(M, mesh->geometry().dofmap(), _x, constants,
79 coefficients);
80 }
81}
82
88template <dolfinx::scalar T, std::floating_point U>
89T assemble_scalar(const Form<T, U>& M)
90{
91 const std::vector<T> constants = pack_constants(M);
92 auto coefficients = allocate_coefficient_storage(M);
93 pack_coefficients(M, coefficients);
94 return assemble_scalar(M, std::span(constants),
95 make_coefficients_span(coefficients));
96}
97
98// -- Vectors ----------------------------------------------------------------
99
110template <dolfinx::scalar T, std::floating_point U>
111void assemble_vector(
112 std::span<T> b, const Form<T, U>& L, std::span<const T> constants,
113 const std::map<std::pair<IntegralType, int>,
114 std::pair<std::span<const T>, int>>& coefficients)
115{
116 impl::assemble_vector(b, L, constants, coefficients);
117}
118
123template <dolfinx::scalar T, std::floating_point U>
124void assemble_vector(std::span<T> b, const Form<T, U>& L)
125{
126 auto coefficients = allocate_coefficient_storage(L);
127 pack_coefficients(L, coefficients);
128 const std::vector<T> constants = pack_constants(L);
129 assemble_vector(b, L, std::span(constants),
130 make_coefficients_span(coefficients));
131}
132
133// FIXME: clarify how x0 is used
134// FIXME: if bcs entries are set
135
136// FIXME: need to pass an array of Vec for x0?
137// FIXME: clarify zeroing of vector
138
151template <dolfinx::scalar T, std::floating_point U>
152void apply_lifting(
153 std::span<T> b, const std::vector<std::shared_ptr<const Form<T, U>>>& a,
154 const std::vector<std::span<const T>>& constants,
155 const std::vector<std::map<std::pair<IntegralType, int>,
156 std::pair<std::span<const T>, int>>>& coeffs,
157 const std::vector<std::vector<std::shared_ptr<const DirichletBC<T, U>>>>&
158 bcs1,
159 const std::vector<std::span<const T>>& x0, T alpha)
160{
161 // If all forms are null, there is nothing to do
162 if (std::ranges::all_of(a, [](auto ptr) { return ptr == nullptr; }))
163 return;
164
165 impl::apply_lifting<T>(b, a, constants, coeffs, bcs1, x0, alpha);
166}
167
180template <dolfinx::scalar T, std::floating_point U>
181void apply_lifting(
182 std::span<T> b, const std::vector<std::shared_ptr<const Form<T, U>>>& a,
183 const std::vector<std::vector<std::shared_ptr<const DirichletBC<T, U>>>>&
184 bcs1,
185 const std::vector<std::span<const T>>& x0, T alpha)
186{
187 std::vector<
188 std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>>
189 coeffs;
190 std::vector<std::vector<T>> constants;
191 for (auto _a : a)
192 {
193 if (_a)
194 {
195 auto coefficients = allocate_coefficient_storage(*_a);
196 pack_coefficients(*_a, coefficients);
197 coeffs.push_back(coefficients);
198 constants.push_back(pack_constants(*_a));
199 }
200 else
201 {
202 coeffs.emplace_back();
203 constants.emplace_back();
204 }
205 }
206
207 std::vector<std::span<const T>> _constants(constants.begin(),
208 constants.end());
209 std::vector<std::map<std::pair<IntegralType, int>,
210 std::pair<std::span<const T>, int>>>
211 _coeffs;
212 std::ranges::transform(coeffs, std::back_inserter(_coeffs),
213 [](auto& c) { return make_coefficients_span(c); });
214
215 apply_lifting(b, a, _constants, _coeffs, bcs1, x0, alpha);
216}
217
218// -- Matrices ---------------------------------------------------------------
219
232template <dolfinx::scalar T, std::floating_point U>
233void assemble_matrix(
234 la::MatSet<T> auto mat_add, const Form<T, U>& a,
235 std::span<const T> constants,
236 const std::map<std::pair<IntegralType, int>,
237 std::pair<std::span<const T>, int>>& coefficients,
238 std::span<const std::int8_t> dof_marker0,
239 std::span<const std::int8_t> dof_marker1)
240
241{
242 std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
243 assert(mesh);
244 if constexpr (std::is_same_v<U, scalar_value_type_t<T>>)
245 {
246 impl::assemble_matrix(mat_add, a, mesh->geometry().dofmap(),
247 mesh->geometry().x(), constants, coefficients,
248 dof_marker0, dof_marker1);
249 }
250 else
251 {
252 auto x = mesh->geometry().x();
253 std::vector<scalar_value_type_t<T>> _x(x.begin(), x.end());
254 impl::assemble_matrix(mat_add, a, mesh->geometry().dofmap(), _x, constants,
255 coefficients, dof_marker0, dof_marker1);
256 }
257}
258
266template <dolfinx::scalar T, std::floating_point U>
267void assemble_matrix(
268 auto mat_add, const Form<T, U>& a, std::span<const T> constants,
269 const std::map<std::pair<IntegralType, int>,
270 std::pair<std::span<const T>, int>>& coefficients,
271 const std::vector<std::shared_ptr<const DirichletBC<T, U>>>& bcs)
272{
273 // Index maps for dof ranges
274 auto map0 = a.function_spaces().at(0)->dofmap()->index_map;
275 auto map1 = a.function_spaces().at(1)->dofmap()->index_map;
276 auto bs0 = a.function_spaces().at(0)->dofmap()->index_map_bs();
277 auto bs1 = a.function_spaces().at(1)->dofmap()->index_map_bs();
278
279 // Build dof markers
280 std::vector<std::int8_t> dof_marker0, dof_marker1;
281 assert(map0);
282 std::int32_t dim0 = bs0 * (map0->size_local() + map0->num_ghosts());
283 assert(map1);
284 std::int32_t dim1 = bs1 * (map1->size_local() + map1->num_ghosts());
285 for (std::size_t k = 0; k < bcs.size(); ++k)
286 {
287 assert(bcs[k]);
288 assert(bcs[k]->function_space());
289 if (a.function_spaces().at(0)->contains(*bcs[k]->function_space()))
290 {
291 dof_marker0.resize(dim0, false);
292 bcs[k]->mark_dofs(dof_marker0);
293 }
294
295 if (a.function_spaces().at(1)->contains(*bcs[k]->function_space()))
296 {
297 dof_marker1.resize(dim1, false);
298 bcs[k]->mark_dofs(dof_marker1);
299 }
300 }
301
302 // Assemble
303 assemble_matrix(mat_add, a, constants, coefficients, dof_marker0,
304 dof_marker1);
305}
306
312template <dolfinx::scalar T, std::floating_point U>
313void assemble_matrix(
314 auto mat_add, const Form<T, U>& a,
315 const std::vector<std::shared_ptr<const DirichletBC<T, U>>>& bcs)
316{
317 // Prepare constants and coefficients
318 const std::vector<T> constants = pack_constants(a);
319 auto coefficients = allocate_coefficient_storage(a);
320 pack_coefficients(a, coefficients);
321
322 // Assemble
323 assemble_matrix(mat_add, a, std::span(constants),
324 make_coefficients_span(coefficients), bcs);
325}
326
337template <dolfinx::scalar T, std::floating_point U>
338void assemble_matrix(auto mat_add, const Form<T, U>& a,
339 std::span<const std::int8_t> dof_marker0,
340 std::span<const std::int8_t> dof_marker1)
341
342{
343 // Prepare constants and coefficients
344 const std::vector<T> constants = pack_constants(a);
345 auto coefficients = allocate_coefficient_storage(a);
346 pack_coefficients(a, coefficients);
347
348 // Assemble
349 assemble_matrix(mat_add, a, std::span(constants),
350 make_coefficients_span(coefficients), dof_marker0,
351 dof_marker1);
352}
353
365template <dolfinx::scalar T>
366void set_diagonal(auto set_fn, std::span<const std::int32_t> rows,
367 T diagonal = 1.0)
368{
369 for (std::size_t i = 0; i < rows.size(); ++i)
370 {
371 std::span diag_span(&diagonal, 1);
372 set_fn(rows.subspan(i, 1), rows.subspan(i, 1), diag_span);
373 }
374}
375
392template <dolfinx::scalar T, std::floating_point U>
394 auto set_fn, const FunctionSpace<U>& V,
395 const std::vector<std::shared_ptr<const DirichletBC<T, U>>>& bcs,
396 T diagonal = 1.0)
397{
398 for (auto& bc : bcs)
399 {
400 assert(bc);
401 if (V.contains(*bc->function_space()))
402 {
403 const auto [dofs, range] = bc->dof_indices();
404 set_diagonal(set_fn, dofs.first(range), diagonal);
405 }
406 }
407}
408} // namespace dolfinx::fem
Definition petsc.h:32
A representation of finite element variational forms.
Definition Form.h:139
This class represents a finite element function space defined by a mesh, a finite element,...
Definition vtk_utils.h:32
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_matrix_impl.h:26
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 utils.h:1008
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:366
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:36
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u into a single array ready for assembly.
Definition utils.h:1348
void pack_coefficients(const Form< T, U > &form, IntegralType integral_type, int id, std::span< T > c, int cstride)
Pack coefficients of a Form for a given integral type and domain id.
Definition utils.h:1062