DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
assemble_vector_impl.h
1// Copyright (C) 2018-2025 Garth N. Wells and Paul T. Kühner
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 "Constant.h"
10#include "DirichletBC.h"
11#include "DofMap.h"
12#include "Form.h"
13#include "traits.h"
14#include "utils.h"
15#include <algorithm>
16#include <basix/mdspan.hpp>
17#include <cstdint>
18#include <dolfinx/common/IndexMap.h>
19#include <dolfinx/mesh/Geometry.h>
20#include <dolfinx/mesh/Mesh.h>
21#include <dolfinx/mesh/Topology.h>
22#include <functional>
23#include <memory>
24#include <optional>
25#include <span>
26#include <vector>
27
28namespace dolfinx::fem
29{
30template <dolfinx::scalar T, std::floating_point U>
31class DirichletBC;
32
33}
34namespace dolfinx::fem::impl
35{
37using mdspan2_t = md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>;
39
79template <int _bs0 = -1, int _bs1 = -1, typename V,
80 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
81 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
82void _lift_bc_cells(
83 V&& b, mdspan2_t x_dofmap,
84 md::mdspan<const scalar_value_t<T>,
85 md::extents<std::size_t, md::dynamic_extent, 3>>
86 x,
87 FEkernel<T> auto kernel, std::span<const std::int32_t> cells,
88 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap0,
89 fem::DofTransformKernel<T> auto P0,
90 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap1,
91 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
92 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
93 std::span<const std::uint32_t> cell_info0,
94 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
95 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha)
96{
97 if (cells.empty())
98 return;
99
100 const auto [dmap0, bs0, cells0] = dofmap0;
101 const auto [dmap1, bs1, cells1] = dofmap1;
102 assert(_bs0 < 0 or _bs0 == bs0);
103 assert(_bs1 < 0 or _bs1 == bs1);
104
105 const int num_rows = bs0 * dmap0.extent(1);
106 const int num_cols = bs1 * dmap1.extent(1);
107
108 // Data structures used in bc application
109 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
110 std::vector<T> Ae(num_rows * num_cols), be(num_rows);
111 assert(cells0.size() == cells.size());
112 assert(cells1.size() == cells.size());
113 for (std::size_t index = 0; index < cells.size(); ++index)
114 {
115 // Cell index in integration domain mesh, test function mesh, and trial
116 // function mesh
117 std::int32_t c = cells[index];
118 std::int32_t c0 = cells0[index];
119 std::int32_t c1 = cells1[index];
120
121 // Get dof maps for cell
122 auto dofs1 = md::submdspan(dmap1, c1, md::full_extent);
123
124 // Check if bc is applied to cell
125 bool has_bc = false;
126 for (std::size_t j = 0; j < dofs1.size(); ++j)
127 {
128 if constexpr (_bs1 > 0)
129 {
130 for (int k = 0; k < _bs1; ++k)
131 {
132 assert(_bs1 * dofs1[j] + k < (int)bc_markers1.size());
133 if (bc_markers1[_bs1 * dofs1[j] + k])
134 {
135 has_bc = true;
136 break;
137 }
138 }
139 }
140 else
141 {
142 for (int k = 0; k < bs1; ++k)
143 {
144 assert(bs1 * dofs1[j] + k < (int)bc_markers1.size());
145 if (bc_markers1[bs1 * dofs1[j] + k])
146 {
147 has_bc = true;
148 break;
149 }
150 }
151 }
152 }
153
154 if (!has_bc)
155 continue;
156
157 // Get cell coordinates/geometry
158 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
159 for (std::size_t i = 0; i < x_dofs.size(); ++i)
160 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
161
162 // Size data structure for assembly
163 auto dofs0 = md::submdspan(dmap0, c0, md::full_extent);
164
165 std::ranges::fill(Ae, 0);
166 kernel(Ae.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
167 nullptr, nullptr, nullptr);
168 P0(Ae, cell_info0, c0, num_cols);
169 P1T(Ae, cell_info1, c1, num_rows);
170
171 // Size data structure for assembly
172 std::ranges::fill(be, 0);
173 for (std::size_t j = 0; j < dofs1.size(); ++j)
174 {
175 if constexpr (_bs1 > 0)
176 {
177 for (int k = 0; k < _bs1; ++k)
178 {
179 const std::int32_t jj = _bs1 * dofs1[j] + k;
180 assert(jj < (int)bc_markers1.size());
181 if (bc_markers1[jj])
182 {
183 const T bc = bc_values1[jj];
184 const T _x0 = x0.empty() ? 0 : x0[jj];
185 // const T _x0 = 0;
186 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
187 for (int m = 0; m < num_rows; ++m)
188 be[m] -= Ae[m * num_cols + _bs1 * j + k] * alpha * (bc - _x0);
189 }
190 }
191 }
192 else
193 {
194 for (int k = 0; k < bs1; ++k)
195 {
196 const std::int32_t jj = bs1 * dofs1[j] + k;
197 assert(jj < (int)bc_markers1.size());
198 if (bc_markers1[jj])
199 {
200 const T bc = bc_values1[jj];
201 const T _x0 = x0.empty() ? 0 : x0[jj];
202 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
203 for (int m = 0; m < num_rows; ++m)
204 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
205 }
206 }
207 }
208 }
209
210 for (std::size_t i = 0; i < dofs0.size(); ++i)
211 {
212 if constexpr (_bs0 > 0)
213 {
214 for (int k = 0; k < _bs0; ++k)
215 b[_bs0 * dofs0[i] + k] += be[_bs0 * i + k];
216 }
217 else
218 {
219 for (int k = 0; k < bs0; ++k)
220 b[bs0 * dofs0[i] + k] += be[bs0 * i + k];
221 }
222 }
223 }
224}
225
263template <typename V,
264 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
265 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
266void _lift_bc_exterior_facets(
267 V&& b, mdspan2_t x_dofmap,
268 md::mdspan<const scalar_value_t<T>,
269 md::extents<std::size_t, md::dynamic_extent, 3>>
270 x,
271 FEkernel<T> auto kernel,
272 md::mdspan<const std::int32_t,
273 md::extents<std::size_t, md::dynamic_extent, 2>>
274 facets,
275 std::tuple<mdspan2_t, int,
276 md::mdspan<const std::int32_t,
277 md::extents<std::size_t, md::dynamic_extent, 2>>>
278 dofmap0,
279 fem::DofTransformKernel<T> auto P0,
280 std::tuple<mdspan2_t, int,
281 md::mdspan<const std::int32_t,
282 md::extents<std::size_t, md::dynamic_extent, 2>>>
283 dofmap1,
284 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
285 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
286 std::span<const std::uint32_t> cell_info0,
287 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
288 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
289 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
290{
291 if (facets.empty())
292 return;
293
294 const auto [dmap0, bs0, facets0] = dofmap0;
295 const auto [dmap1, bs1, facets1] = dofmap1;
296
297 const int num_rows = bs0 * dmap0.extent(1);
298 const int num_cols = bs1 * dmap1.extent(1);
299
300 // Data structures used in bc application
301 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
302 std::vector<T> Ae(num_rows * num_cols), be(num_rows);
303 assert(facets0.size() == facets.size());
304 assert(facets1.size() == facets.size());
305 for (std::size_t index = 0; index < facets.extent(0); ++index)
306 {
307 // Cell in integration domain, test function and trial function
308 // meshes
309 std::int32_t cell = facets(index, 0);
310 std::int32_t cell0 = facets0(index, 0);
311 std::int32_t cell1 = facets1(index, 0);
312
313 // Local facet index
314 std::int32_t local_facet = facets(index, 1);
315
316 // Get dof maps for cell
317 auto dofs1 = md::submdspan(dmap1, cell1, md::full_extent);
318
319 // Check if bc is applied to cell
320 bool has_bc = false;
321 for (std::size_t j = 0; j < dofs1.size(); ++j)
322 {
323 for (int k = 0; k < bs1; ++k)
324 {
325 if (bc_markers1[bs1 * dofs1[j] + k])
326 {
327 has_bc = true;
328 break;
329 }
330 }
331 }
332
333 if (!has_bc)
334 continue;
335
336 // Get cell coordinates/geometry
337 auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
338 for (std::size_t i = 0; i < x_dofs.size(); ++i)
339 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
340
341 // Size data structure for assembly
342 auto dofs0 = md::submdspan(dmap0, cell0, md::full_extent);
343
344 // Permutations
345 std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
346 std::ranges::fill(Ae, 0);
347 kernel(Ae.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
348 &local_facet, &perm, nullptr);
349 P0(Ae, cell_info0, cell0, num_cols);
350 P1T(Ae, cell_info1, cell1, num_rows);
351
352 // Size data structure for assembly
353 std::ranges::fill(be, 0);
354 for (std::size_t j = 0; j < dofs1.size(); ++j)
355 {
356 for (int k = 0; k < bs1; ++k)
357 {
358 const std::int32_t jj = bs1 * dofs1[j] + k;
359 if (bc_markers1[jj])
360 {
361 const T bc = bc_values1[jj];
362 const T _x0 = x0.empty() ? 0 : x0[jj];
363 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
364 for (int m = 0; m < num_rows; ++m)
365 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
366 }
367 }
368 }
369
370 for (std::size_t i = 0; i < dofs0.size(); ++i)
371 for (int k = 0; k < bs0; ++k)
372 b[bs0 * dofs0[i] + k] += be[bs0 * i + k];
373 }
374}
375
418template <typename V,
419 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
420 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
421void _lift_bc_interior_facets(
422 V&& b, mdspan2_t x_dofmap,
423 md::mdspan<const scalar_value_t<T>,
424 md::extents<std::size_t, md::dynamic_extent, 3>>
425 x,
426 FEkernel<T> auto kernel,
427 md::mdspan<const std::int32_t,
428 md::extents<std::size_t, md::dynamic_extent, 2, 2>>
429 facets,
430 std::tuple<mdspan2_t, int,
431 md::mdspan<const std::int32_t,
432 md::extents<std::size_t, md::dynamic_extent, 2, 2>>>
433 dofmap0,
434 fem::DofTransformKernel<T> auto P0,
435 std::tuple<mdspan2_t, int,
436 md::mdspan<const std::int32_t,
437 md::extents<std::size_t, md::dynamic_extent, 2, 2>>>
438 dofmap1,
439 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
440 md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
441 md::dynamic_extent>>
442 coeffs,
443 std::span<const std::uint32_t> cell_info0,
444 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
445 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
446 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
447{
448 if (facets.empty())
449 return;
450
451 const auto [dmap0, bs0, facets0] = dofmap0;
452 const auto [dmap1, bs1, facets1] = dofmap1;
453
454 const int num_dofs0 = dmap0.extent(1);
455 const int num_dofs1 = dmap1.extent(1);
456 const int num_rows = bs0 * 2 * num_dofs0;
457 const int num_cols = bs1 * 2 * num_dofs1;
458
459 // Data structures used in assembly
460 using X = scalar_value_t<T>;
461 std::vector<X> cdofs(2 * x_dofmap.extent(1) * 3);
462 std::span<X> cdofs0(cdofs.data(), x_dofmap.extent(1) * 3);
463 std::span<X> cdofs1(cdofs.data() + x_dofmap.extent(1) * 3,
464 x_dofmap.extent(1) * 3);
465 std::vector<T> Ae(num_rows * num_cols), be(num_rows);
466
467 // Temporaries for joint dofmaps
468 std::vector<std::int32_t> dmapjoint0(2 * num_dofs0);
469 std::vector<std::int32_t> dmapjoint1(2 * num_dofs1);
470
471 assert(facets0.size() == facets.size());
472 assert(facets1.size() == facets.size());
473 for (std::size_t f = 0; f < facets.extent(0); ++f)
474 {
475 // Cells in integration domain, test function domain and trial
476 // function domain meshes
477 std::array cells{facets(f, 0, 0), facets(f, 1, 0)};
478 std::array cells0{facets0(f, 0, 0), facets0(f, 1, 0)};
479 std::array cells1{facets1(f, 0, 0), facets1(f, 1, 0)};
480
481 // Local facet indices
482 std::array local_facet = {facets(f, 0, 1), facets(f, 1, 1)};
483
484 // Get cell geometry
485 auto x_dofs0 = md::submdspan(x_dofmap, cells[0], md::full_extent);
486 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
487 std::copy_n(&x(x_dofs0[i], 0), 3, std::next(cdofs0.begin(), 3 * i));
488 auto x_dofs1 = md::submdspan(x_dofmap, cells[1], md::full_extent);
489 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
490 std::copy_n(&x(x_dofs1[i], 0), 3, std::next(cdofs1.begin(), 3 * i));
491
492 // Get dof maps for cells and pack
493 // When integrating over interfaces between two domains, the test function
494 // might only be defined on one side, so we check which cells exist in the
495 // test function domain
496 std::span dmap0_cell0
497 = cells0[0] >= 0
498 ? std::span(dmap0.data_handle() + cells0[0] * num_dofs0,
499 num_dofs0)
500 : std::span<const std::int32_t>();
501 std::span dmap0_cell1
502 = cells0[1] >= 0
503 ? std::span(dmap0.data_handle() + cells0[1] * num_dofs0,
504 num_dofs0)
505 : std::span<const std::int32_t>();
506
507 std::ranges::copy(dmap0_cell0, dmapjoint0.begin());
508 std::ranges::copy(dmap0_cell1, std::next(dmapjoint0.begin(), num_dofs0));
509
510 // Check which cells exist in the trial function domain
511 std::span<const std::int32_t> dmap1_cell0
512 = cells1[0] >= 0
513 ? std::span(dmap1.data_handle() + cells1[0] * num_dofs1,
514 num_dofs1)
515 : std::span<const std::int32_t>();
516 std::span<const std::int32_t> dmap1_cell1
517 = cells1[1] >= 0
518 ? std::span(dmap1.data_handle() + cells1[1] * num_dofs1,
519 num_dofs1)
520 : std::span<const std::int32_t>();
521
522 std::ranges::copy(dmap1_cell0, dmapjoint1.begin());
523 std::ranges::copy(dmap1_cell1, std::next(dmapjoint1.begin(), num_dofs1));
524
525 // Check if bc is applied to cell0
526 bool has_bc = false;
527 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
528 {
529 for (int k = 0; k < bs1; ++k)
530 {
531 if (bc_markers1[bs1 * dmap1_cell0[j] + k])
532 {
533 has_bc = true;
534 break;
535 }
536 }
537 }
538
539 // Check if bc is applied to cell1
540 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
541 {
542 for (int k = 0; k < bs1; ++k)
543 {
544 if (bc_markers1[bs1 * dmap1_cell1[j] + k])
545 {
546 has_bc = true;
547 break;
548 }
549 }
550 }
551
552 if (!has_bc)
553 continue;
554
555 // Tabulate tensor
556 std::ranges::fill(Ae, 0);
557 std::array perm = perms.empty()
558 ? std::array<std::uint8_t, 2>{0, 0}
559 : std::array{perms(cells[0], local_facet[0]),
560 perms(cells[1], local_facet[1])};
561 kernel(Ae.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
562 local_facet.data(), perm.data(), nullptr);
563
564 if (cells0[0] >= 0)
565 P0(Ae, cell_info0, cells0[0], num_cols);
566 if (cells0[1] >= 0)
567 {
568 std::span sub_Ae0(Ae.data() + bs0 * num_dofs0 * num_cols,
569 bs0 * num_dofs1 * num_cols);
570 P0(sub_Ae0, cell_info0, cells0[1], num_cols);
571 }
572 if (cells1[0] >= 0)
573 P1T(Ae, cell_info1, cells1[0], num_rows);
574
575 if (cells1[1] >= 0)
576 {
577 for (int row = 0; row < num_rows; ++row)
578 {
579 // DOFs for dmap1 and cell1 are not stored contiguously in
580 // the block matrix, so each row needs a separate span access
581 std::span sub_Ae1(Ae.data() + row * num_cols + bs1 * num_dofs1,
582 bs1 * num_dofs1);
583 P1T(sub_Ae1, cell_info1, cells1[1], 1);
584 }
585 }
586
587 std::ranges::fill(be, 0);
588
589 // Compute b = b - A*b for cell0
590 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
591 {
592 for (int k = 0; k < bs1; ++k)
593 {
594 const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
595 if (bc_markers1[jj])
596 {
597 const T bc = bc_values1[jj];
598 const T _x0 = x0.empty() ? 0 : x0[jj];
599 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
600 for (int m = 0; m < num_rows; ++m)
601 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
602 }
603 }
604 }
605
606 // Compute b = b - A*b for cell1
607 const int offset = bs1 * num_dofs1;
608 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
609 {
610 for (int k = 0; k < bs1; ++k)
611 {
612 const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
613 if (bc_markers1[jj])
614 {
615 const T bc = bc_values1[jj];
616 const T _x0 = x0.empty() ? 0 : x0[jj];
617 // be -= Ae.col(offset + bs1 * j + k) * alpha * (bc - x0[jj]);
618 for (int m = 0; m < num_rows; ++m)
619 {
620 be[m]
621 -= Ae[m * num_cols + offset + bs1 * j + k] * alpha * (bc - _x0);
622 }
623 }
624 }
625 }
626
627 for (std::size_t i = 0; i < dmap0_cell0.size(); ++i)
628 for (int k = 0; k < bs0; ++k)
629 b[bs0 * dmap0_cell0[i] + k] += be[bs0 * i + k];
630
631 const int offset_be = bs0 * num_dofs0;
632 for (std::size_t i = 0; i < dmap0_cell1.size(); ++i)
633 for (int k = 0; k < bs0; ++k)
634 b[bs0 * dmap0_cell1[i] + k] += be[offset_be + bs0 * i + k];
635 }
636}
637
661template <int _bs = -1, typename V,
662 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
663 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
664
665void assemble_cells(
666 fem::DofTransformKernel<T> auto P0, V&& b, mdspan2_t x_dofmap,
667 md::mdspan<const scalar_value_t<T>,
668 md::extents<std::size_t, md::dynamic_extent, 3>>
669 x,
670 std::span<const std::int32_t> cells,
671 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap,
672 FEkernel<T> auto kernel, std::span<const T> constants,
673 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
674 std::span<const std::uint32_t> cell_info0)
675{
676 if (cells.empty())
677 return;
678
679 const auto [dmap, bs, cells0] = dofmap;
680 assert(_bs < 0 or _bs == bs);
681
682 // Create data structures used in assembly
683 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
684 std::vector<T> be(bs * dmap.extent(1));
685
686 // Iterate over active cells
687 for (std::size_t index = 0; index < cells.size(); ++index)
688 {
689 // Integration domain celland test function cell
690 std::int32_t c = cells[index];
691 std::int32_t c0 = cells0[index];
692
693 // Get cell coordinates/geometry
694 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
695 for (std::size_t i = 0; i < x_dofs.size(); ++i)
696 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
697
698 // Tabulate vector for cell
699 std::ranges::fill(be, 0);
700 kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
701 nullptr, nullptr, nullptr);
702 P0(be, cell_info0, c0, 1);
703
704 // Scatter cell vector to 'global' vector array
705 auto dofs = md::submdspan(dmap, c0, md::full_extent);
706 if constexpr (_bs > 0)
707 {
708 for (std::size_t i = 0; i < dofs.size(); ++i)
709 for (int k = 0; k < _bs; ++k)
710 b[_bs * dofs[i] + k] += be[_bs * i + k];
711 }
712 else
713 {
714 for (std::size_t i = 0; i < dofs.size(); ++i)
715 for (int k = 0; k < bs; ++k)
716 b[bs * dofs[i] + k] += be[bs * i + k];
717 }
718 }
719}
720
745template <int _bs = -1, typename V,
746 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
747 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
748void assemble_exterior_facets(
749 fem::DofTransformKernel<T> auto P0, V&& b, mdspan2_t x_dofmap,
750 md::mdspan<const scalar_value_t<T>,
751 md::extents<std::size_t, md::dynamic_extent, 3>>
752 x,
753 md::mdspan<const std::int32_t,
754 std::extents<std::size_t, md::dynamic_extent, 2>>
755 facets,
756 std::tuple<mdspan2_t, int,
757 md::mdspan<const std::int32_t,
758 std::extents<std::size_t, md::dynamic_extent, 2>>>
759 dofmap,
760 FEkernel<T> auto kernel, std::span<const T> constants,
761 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
762 std::span<const std::uint32_t> cell_info0,
763 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
764{
765 if (facets.empty())
766 return;
767
768 const auto [dmap, bs, facets0] = dofmap;
769 assert(_bs < 0 or _bs == bs);
770
771 // Create data structures used in assembly
772 const int num_dofs = dmap.extent(1);
773 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
774 std::vector<T> be(bs * num_dofs);
775 assert(facets0.size() == facets.size());
776 for (std::size_t f = 0; f < facets.extent(0); ++f)
777 {
778 // Cell in the integration domain, local facet index relative to the
779 // integration domain cell, and cell in the test function mesh
780 std::int32_t cell = facets(f, 0);
781 std::int32_t local_facet = facets(f, 1);
782 std::int32_t cell0 = facets0(f, 0);
783
784 // Get cell coordinates/geometry
785 auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
786 for (std::size_t i = 0; i < x_dofs.size(); ++i)
787 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
788
789 // Permutations
790 std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
791
792 // Tabulate element vector
793 std::ranges::fill(be, 0);
794 kernel(be.data(), &coeffs(f, 0), constants.data(), cdofs.data(),
795 &local_facet, &perm, nullptr);
796 P0(be, cell_info0, cell0, 1);
797
798 // Add element vector to global vector
799 auto dofs = md::submdspan(dmap, cell0, md::full_extent);
800 if constexpr (_bs > 0)
801 {
802 for (std::size_t i = 0; i < dofs.size(); ++i)
803 for (int k = 0; k < _bs; ++k)
804 b[_bs * dofs[i] + k] += be[_bs * i + k];
805 }
806 else
807 {
808 for (std::size_t i = 0; i < dofs.size(); ++i)
809 for (int k = 0; k < bs; ++k)
810 b[bs * dofs[i] + k] += be[bs * i + k];
811 }
812 }
813}
814
840template <int _bs = -1, typename V,
841 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
842 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
843void assemble_interior_facets(
844 fem::DofTransformKernel<T> auto P0, V&& b, mdspan2_t x_dofmap,
845 md::mdspan<const scalar_value_t<T>,
846 md::extents<std::size_t, md::dynamic_extent, 3>>
847 x,
848 md::mdspan<const std::int32_t,
849 std::extents<std::size_t, md::dynamic_extent, 2, 2>>
850 facets,
851 std::tuple<const DofMap&, int,
852 md::mdspan<const std::int32_t,
853 std::extents<std::size_t, md::dynamic_extent, 2, 2>>>
854 dofmap,
855 FEkernel<T> auto kernel, std::span<const T> constants,
856 md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
857 md::dynamic_extent>>
858 coeffs,
859 std::span<const std::uint32_t> cell_info0,
860 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
861{
862 using X = scalar_value_t<T>;
863
864 if (facets.empty())
865 return;
866
867 const auto [dmap, bs, facets0] = dofmap;
868 assert(_bs < 0 or _bs == bs);
869
870 // Create data structures used in assembly
871 std::vector<X> cdofs(2 * x_dofmap.extent(1) * 3);
872 std::span<X> cdofs0(cdofs.data(), x_dofmap.extent(1) * 3);
873 std::span<X> cdofs1(cdofs.data() + x_dofmap.extent(1) * 3,
874 x_dofmap.extent(1) * 3);
875
876 const std::size_t dmap_size = dmap.map().extent(1);
877 std::vector<T> be(bs * 2 * dmap_size);
878
879 assert(facets0.size() == facets.size());
880 for (std::size_t f = 0; f < facets.extent(0); ++f)
881 {
882 // Cells in integration domain and test function domain meshes
883 std::array<std::int32_t, 2> cells{facets(f, 0, 0), facets(f, 1, 0)};
884 std::array<std::int32_t, 2> cells0{facets0(f, 0, 0), facets0(f, 1, 0)};
885
886 // Local facet indices
887 std::array<std::int32_t, 2> local_facet{facets(f, 0, 1), facets(f, 1, 1)};
888
889 // Get cell geometry
890 auto x_dofs0 = md::submdspan(x_dofmap, cells[0], md::full_extent);
891 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
892 std::copy_n(&x(x_dofs0[i], 0), 3, std::next(cdofs0.begin(), 3 * i));
893 auto x_dofs1 = md::submdspan(x_dofmap, cells[1], md::full_extent);
894 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
895 std::copy_n(&x(x_dofs1[i], 0), 3, std::next(cdofs1.begin(), 3 * i));
896
897 // Get dofmaps for cells. When integrating over interfaces between
898 // two domains, the test function might only be defined on one side,
899 // so we check which cells exist in the test function domain.
900 std::span dmap0 = cells0[0] >= 0 ? dmap.cell_dofs(cells0[0])
901 : std::span<const std::int32_t>();
902 std::span dmap1 = cells0[1] >= 0 ? dmap.cell_dofs(cells0[1])
903 : std::span<const std::int32_t>();
904
905 // Tabulate element vector
906 std::ranges::fill(be, 0);
907 std::array perm = perms.empty()
908 ? std::array<std::uint8_t, 2>{0, 0}
909 : std::array{perms(cells[0], local_facet[0]),
910 perms(cells[1], local_facet[1])};
911 kernel(be.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
912 local_facet.data(), perm.data(), nullptr);
913
914 if (cells0[0] >= 0)
915 P0(be, cell_info0, cells0[0], 1);
916 if (cells0[1] >= 0)
917 {
918 std::span sub_be(be.data() + bs * dmap_size, bs * dmap_size);
919 P0(sub_be, cell_info0, cells0[1], 1);
920 }
921
922 // Add element vector to global vector
923 if constexpr (_bs > 0)
924 {
925 for (std::size_t i = 0; i < dmap0.size(); ++i)
926 for (int k = 0; k < _bs; ++k)
927 b[_bs * dmap0[i] + k] += be[_bs * i + k];
928 for (std::size_t i = 0; i < dmap1.size(); ++i)
929 for (int k = 0; k < _bs; ++k)
930 b[_bs * dmap1[i] + k] += be[_bs * (i + dmap_size) + k];
931 }
932 else
933 {
934 for (std::size_t i = 0; i < dmap0.size(); ++i)
935 for (int k = 0; k < bs; ++k)
936 b[bs * dmap0[i] + k] += be[bs * i + k];
937 for (std::size_t i = 0; i < dmap1.size(); ++i)
938 for (int k = 0; k < bs; ++k)
939 b[bs * dmap1[i] + k] += be[bs * (i + dmap_size) + k];
940 }
941 }
942}
943
969template <dolfinx::scalar T, int _bs = -1>
970void assemble_vertices(
971 fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
972 md::mdspan<const scalar_value_t<T>,
973 md::extents<std::size_t, md::dynamic_extent, 3>>
974 x,
975 md::mdspan<const std::int32_t,
976 md::extents<std::size_t, md::dynamic_extent, 2>>
977 vertices,
978 std::tuple<mdspan2_t, int,
979 md::mdspan<const std::int32_t,
980 md::extents<std::size_t, md::dynamic_extent, 2>>>
981 dofmap,
982 FEkernel<T> auto kernel, std::span<const T> constants,
983 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
984 std::span<const std::uint32_t> cell_info0)
985{
986 if (vertices.empty())
987 return;
988
989 const auto [dmap, bs, vertices0] = dofmap;
990 assert(_bs < 0 or _bs == bs);
991
992 // Create data structures used in assembly
993 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
994 std::vector<T> be(bs * dmap.extent(1));
995
996 // Iterate over active vertices
997 for (std::size_t index = 0; index < vertices.extent(0); ++index)
998 {
999 // Integration domain cell, local index, and test function cell
1000 std::int32_t cell = vertices(index, 0);
1001 std::int32_t local_index = vertices(index, 1);
1002 std::int32_t c0 = vertices0(index, 0);
1003
1004 // Get cell coordinates/geometry
1005 auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
1006 for (std::size_t i = 0; i < x_dofs.size(); ++i)
1007 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
1008
1009 // Tabulate vector for vertex
1010 std::ranges::fill(be, 0);
1011 kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
1012 &local_index, nullptr, nullptr);
1013 P0(be, cell_info0, c0, 1);
1014
1015 // Scatter vertex vector to 'global' vector array
1016 auto dofs = md::submdspan(dmap, c0, md::full_extent);
1017 if constexpr (_bs > 0)
1018 {
1019 for (std::size_t i = 0; i < dofs.size(); ++i)
1020 for (int k = 0; k < _bs; ++k)
1021 b[_bs * dofs[i] + k] += be[_bs * i + k];
1022 }
1023 else
1024 {
1025 for (std::size_t i = 0; i < dofs.size(); ++i)
1026 for (int k = 0; k < bs; ++k)
1027 b[bs * dofs[i] + k] += be[bs * i + k];
1028 }
1029 }
1030}
1031
1048template <typename V, std::floating_point U,
1049 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
1050 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
1051void lift_bc(V&& b, const Form<T, U>& a, mdspan2_t x_dofmap,
1052 md::mdspan<const scalar_value_t<T>,
1053 md::extents<std::size_t, md::dynamic_extent, 3>>
1054 x,
1055 std::span<const T> constants,
1056 const std::map<std::pair<IntegralType, int>,
1057 std::pair<std::span<const T>, int>>& coefficients,
1058 std::span<const T> bc_values1,
1059 std::span<const std::int8_t> bc_markers1, std::span<const T> x0,
1060 T alpha)
1061{
1062 // Integration domain mesh
1063 std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
1064 assert(mesh);
1065
1066 // Test function mesh
1067 auto mesh0 = a.function_spaces().at(0)->mesh();
1068 assert(mesh0);
1069
1070 // Trial function mesh
1071 auto mesh1 = a.function_spaces().at(1)->mesh();
1072 assert(mesh1);
1073
1074 // Get dofmap for columns and rows of a
1075 assert(a.function_spaces().at(0));
1076 assert(a.function_spaces().at(1));
1077 auto dofmap0 = a.function_spaces()[0]->dofmap()->map();
1078 const int bs0 = a.function_spaces()[0]->dofmap()->bs();
1079 auto element0 = a.function_spaces()[0]->element();
1080 auto dofmap1 = a.function_spaces()[1]->dofmap()->map();
1081 const int bs1 = a.function_spaces()[1]->dofmap()->bs();
1082 auto element1 = a.function_spaces()[1]->element();
1083 assert(element0);
1084
1085 std::span<const std::uint32_t> cell_info0;
1086 std::span<const std::uint32_t> cell_info1;
1087 // TODO: Check for each element instead
1088 if (element0->needs_dof_transformations()
1089 or element1->needs_dof_transformations() or a.needs_facet_permutations())
1090 {
1091 mesh0->topology_mutable()->create_entity_permutations();
1092 mesh1->topology_mutable()->create_entity_permutations();
1093 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
1094 cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
1095 }
1096
1097 fem::DofTransformKernel<T> auto P0
1098 = element0->template dof_transformation_fn<T>(doftransform::standard);
1099 fem::DofTransformKernel<T> auto P1T
1100 = element1->template dof_transformation_right_fn<T>(
1102
1103 for (int i = 0; i < a.num_integrals(IntegralType::cell, 0); ++i)
1104 {
1105 auto kernel = a.kernel(IntegralType::cell, i, 0);
1106 assert(kernel);
1107 auto& [_coeffs, cstride] = coefficients.at({IntegralType::cell, i});
1108 std::span cells = a.domain(IntegralType::cell, i, 0);
1109 std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, 0);
1110 std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, 0);
1111 assert(_coeffs.size() == cells.size() * cstride);
1112 auto coeffs = md::mdspan(_coeffs.data(), cells.size(), cstride);
1113 if (bs0 == 1 and bs1 == 1)
1114 {
1115 _lift_bc_cells<1, 1>(b, x_dofmap, x, kernel, cells,
1116 {dofmap0, bs0, cells0}, P0, {dofmap1, bs1, cells1},
1117 P1T, constants, coeffs, cell_info0, cell_info1,
1118 bc_values1, bc_markers1, x0, alpha);
1119 }
1120 else if (bs0 == 3 and bs1 == 3)
1121 {
1122 _lift_bc_cells<3, 3>(b, x_dofmap, x, kernel, cells,
1123 {dofmap0, bs0, cells0}, P0, {dofmap1, bs1, cells1},
1124 P1T, constants, coeffs, cell_info0, cell_info1,
1125 bc_values1, bc_markers1, x0, alpha);
1126 }
1127 else
1128 {
1129 _lift_bc_cells(b, x_dofmap, x, kernel, cells, {dofmap0, bs0, cells0}, P0,
1130 {dofmap1, bs1, cells1}, P1T, constants, coeffs, cell_info0,
1131 cell_info1, bc_values1, bc_markers1, x0, alpha);
1132 }
1133 }
1134
1135 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
1136 if (a.needs_facet_permutations())
1137 {
1138 mesh::CellType cell_type = mesh->topology()->cell_type();
1139 int num_facets_per_cell
1140 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
1141 mesh->topology_mutable()->create_entity_permutations();
1142 const std::vector<std::uint8_t>& p
1143 = mesh->topology()->get_facet_permutations();
1144 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
1145 num_facets_per_cell);
1146 }
1147
1148 for (int i = 0; i < a.num_integrals(IntegralType::exterior_facet, 0); ++i)
1149 {
1150 auto kernel = a.kernel(IntegralType::exterior_facet, i, 0);
1151 assert(kernel);
1152 auto& [coeffs, cstride]
1153 = coefficients.at({IntegralType::exterior_facet, i});
1154
1155 using mdspanx2_t
1156 = md::mdspan<const std::int32_t,
1157 md::extents<std::size_t, md::dynamic_extent, 2>>;
1158 std::span f = a.domain(IntegralType::exterior_facet, i, 0);
1159 mdspanx2_t facets(f.data(), f.size() / 2, 2);
1160 std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0);
1161 mdspanx2_t facets0(f0.data(), f0.size() / 2, 2);
1162 std::span f1 = a.domain_arg(IntegralType::exterior_facet, 1, i, 0);
1163 mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
1164 assert(coeffs.size() == facets.extent(0) * cstride);
1165 _lift_bc_exterior_facets(
1166 b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
1167 {dofmap1, bs1, facets1}, P1T, constants,
1168 md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
1169 cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
1170 }
1171
1172 for (int i = 0; i < a.num_integrals(IntegralType::interior_facet, 0); ++i)
1173 {
1174 auto kernel = a.kernel(IntegralType::interior_facet, i, 0);
1175 assert(kernel);
1176 auto& [coeffs, cstride]
1177 = coefficients.at({IntegralType::interior_facet, i});
1178
1179 using mdspanx22_t
1180 = md::mdspan<const std::int32_t,
1181 md::extents<std::size_t, md::dynamic_extent, 2, 2>>;
1182 using mdspanx2x_t
1183 = md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
1184 md::dynamic_extent>>;
1185 std::span f = a.domain(IntegralType::interior_facet, i, 0);
1186 mdspanx22_t facets(f.data(), f.size() / 4, 2, 2);
1187 std::span f0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0);
1188 mdspanx22_t facets0(f0.data(), f0.size() / 4, 2, 2);
1189 std::span f1 = a.domain_arg(IntegralType::interior_facet, 1, i, 0);
1190 mdspanx22_t facets1(f1.data(), f1.size() / 4, 2, 2);
1191 _lift_bc_interior_facets(
1192 b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
1193 {dofmap1, bs1, facets1}, P1T, constants,
1194 mdspanx2x_t(coeffs.data(), facets.extent(0), 2, cstride), cell_info0,
1195 cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
1196 }
1197}
1198
1219template <typename V, std::floating_point U,
1220 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
1221 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
1222void apply_lifting(
1223 V&& b,
1224 std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> a,
1225 const std::vector<std::span<const T>>& constants,
1226 const std::vector<std::map<std::pair<IntegralType, int>,
1227 std::pair<std::span<const T>, int>>>& coeffs,
1228 const std::vector<
1229 std::vector<std::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
1230 const std::vector<std::span<const T>>& x0, T alpha)
1231{
1232 if (!x0.empty() and x0.size() != a.size())
1233 {
1234 throw std::runtime_error(
1235 "Mismatch in size between x0 and bilinear form in assembler.");
1236 }
1237
1238 if (a.size() != bcs1.size())
1239 {
1240 throw std::runtime_error(
1241 "Mismatch in size between a and bcs in assembler.");
1242 }
1243
1244 for (std::size_t j = 0; j < a.size(); ++j)
1245 {
1246 std::vector<std::int8_t> bc_markers1;
1247 std::vector<T> bc_values1;
1248 if (a[j] and !bcs1[j].empty())
1249 {
1250 // Extract data from mesh
1251 std::shared_ptr<const mesh::Mesh<U>> mesh = a[j]->get().mesh();
1252 if (!mesh)
1253 throw std::runtime_error("Unable to extract a mesh.");
1254 mdspan2_t x_dofmap = mesh->geometry().dofmap();
1255 std::span _x = mesh->geometry().x();
1256 md::mdspan<const scalar_value_t<T>,
1257 md::extents<std::size_t, md::dynamic_extent, 3>>
1258 x(_x.data(), _x.size() / 3, 3);
1259
1260 assert(a[j]->get().function_spaces().at(0));
1261 auto V1 = a[j]->get().function_spaces()[1];
1262 assert(V1);
1263 auto map1 = V1->dofmap()->index_map;
1264 const int bs1 = V1->dofmap()->index_map_bs();
1265 assert(map1);
1266 const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
1267 bc_markers1.assign(crange, false);
1268 bc_values1.assign(crange, 0);
1269 for (auto& bc : bcs1[j])
1270 {
1271 bc.get().mark_dofs(bc_markers1);
1272 bc.get().set(bc_values1, std::nullopt, 1);
1273 }
1274
1275 if (!x0.empty())
1276 {
1277 lift_bc(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j],
1278 std::span<const T>(bc_values1), bc_markers1, x0[j], alpha);
1279 }
1280 else
1281 {
1282 lift_bc(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j],
1283 std::span<const T>(bc_values1), bc_markers1,
1284 std::span<const T>(), alpha);
1285 }
1286 }
1287 }
1288}
1289
1297template <typename V, std::floating_point U,
1298 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
1299 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
1300void assemble_vector(
1301 V&& b, const Form<T, U>& L,
1302 md::mdspan<const scalar_value_t<T>,
1303 md::extents<std::size_t, md::dynamic_extent, 3>>
1304 x,
1305 std::span<const T> constants,
1306 const std::map<std::pair<IntegralType, int>,
1307 std::pair<std::span<const T>, int>>& coefficients)
1308{
1309 // Integration domain mesh
1310 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1311 assert(mesh);
1312
1313 // Test function mesh
1314 auto mesh0 = L.function_spaces().at(0)->mesh();
1315 assert(mesh0);
1316
1317 const int num_cell_types = mesh->topology()->cell_types().size();
1318 for (int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
1319 {
1320 // Geometry dofmap and data
1321 mdspan2_t x_dofmap = mesh->geometry().dofmap(cell_type_idx);
1322
1323 // Get dofmap data
1324 assert(L.function_spaces().at(0));
1325 auto element = L.function_spaces().at(0)->elements(cell_type_idx);
1326 assert(element);
1327 std::shared_ptr<const fem::DofMap> dofmap
1328 = L.function_spaces().at(0)->dofmaps(cell_type_idx);
1329 assert(dofmap);
1330 auto dofs = dofmap->map();
1331 const int bs = dofmap->bs();
1332
1333 fem::DofTransformKernel<T> auto P0
1334 = element->template dof_transformation_fn<T>(doftransform::standard);
1335
1336 std::span<const std::uint32_t> cell_info0;
1337 if (element->needs_dof_transformations() or L.needs_facet_permutations())
1338 {
1339 mesh0->topology_mutable()->create_entity_permutations();
1340 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
1341 }
1342
1343 for (int i = 0; i < L.num_integrals(IntegralType::cell, 0); ++i)
1344 {
1345 auto fn = L.kernel(IntegralType::cell, i, cell_type_idx);
1346 assert(fn);
1347 std::span cells = L.domain(IntegralType::cell, i, cell_type_idx);
1348 std::span cells0 = L.domain_arg(IntegralType::cell, 0, i, cell_type_idx);
1349 auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
1350 assert(cells.size() * cstride == coeffs.size());
1351 if (bs == 1)
1352 {
1353 impl::assemble_cells<1>(
1354 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1355 md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
1356 }
1357 else if (bs == 3)
1358 {
1359 impl::assemble_cells<3>(
1360 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1361 md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
1362 }
1363 else
1364 {
1365 impl::assemble_cells(
1366 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1367 md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
1368 }
1369 }
1370
1371 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
1372 if (L.needs_facet_permutations())
1373 {
1374 mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx];
1375 int num_facets_per_cell
1376 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
1377 mesh->topology_mutable()->create_entity_permutations();
1378 const std::vector<std::uint8_t>& p
1379 = mesh->topology()->get_facet_permutations();
1380 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
1381 num_facets_per_cell);
1382 }
1383
1384 using mdspanx2_t
1385 = md::mdspan<const std::int32_t,
1386 md::extents<std::size_t, md::dynamic_extent, 2>>;
1387
1388 for (int i = 0; i < L.num_integrals(IntegralType::exterior_facet, 0); ++i)
1389 {
1390 auto fn = L.kernel(IntegralType::exterior_facet, i, 0);
1391 assert(fn);
1392 auto& [coeffs, cstride]
1393 = coefficients.at({IntegralType::exterior_facet, i});
1394 std::span f = L.domain(IntegralType::exterior_facet, i, 0);
1395 mdspanx2_t facets(f.data(), f.size() / 2, 2);
1396 std::span f1 = L.domain_arg(IntegralType::exterior_facet, 0, i, 0);
1397 mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
1398 assert((facets.size() / 2) * cstride == coeffs.size());
1399 if (bs == 1)
1400 {
1401 impl::assemble_exterior_facets<1>(
1402 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1403 md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
1404 perms);
1405 }
1406 else if (bs == 3)
1407 {
1408 impl::assemble_exterior_facets<3>(
1409 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1410 md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
1411 perms);
1412 }
1413 else
1414 {
1415 impl::assemble_exterior_facets(
1416 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1417 md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
1418 perms);
1419 }
1420 }
1421
1422 for (int i = 0; i < L.num_integrals(IntegralType::interior_facet, 0); ++i)
1423 {
1424 using mdspanx22_t
1425 = md::mdspan<const std::int32_t,
1426 md::extents<std::size_t, md::dynamic_extent, 2, 2>>;
1427 using mdspanx2x_t
1428 = md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
1429 md::dynamic_extent>>;
1430
1431 auto fn = L.kernel(IntegralType::interior_facet, i, 0);
1432 assert(fn);
1433 auto& [coeffs, cstride]
1434 = coefficients.at({IntegralType::interior_facet, i});
1435 std::span facets = L.domain(IntegralType::interior_facet, i, 0);
1436 std::span facets1 = L.domain_arg(IntegralType::interior_facet, 0, i, 0);
1437 assert((facets.size() / 4) * 2 * cstride == coeffs.size());
1438 if (bs == 1)
1439 {
1440 impl::assemble_interior_facets<1>(
1441 P0, b, x_dofmap, x,
1442 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1443 {*dofmap, bs,
1444 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1445 fn, constants,
1446 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1447 cell_info0, perms);
1448 }
1449 else if (bs == 3)
1450 {
1451 impl::assemble_interior_facets<3>(
1452 P0, b, x_dofmap, x,
1453 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1454 {*dofmap, bs,
1455 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1456 fn, constants,
1457 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1458 cell_info0, perms);
1459 }
1460 else
1461 {
1462 impl::assemble_interior_facets(
1463 P0, b, x_dofmap, x,
1464 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1465 {*dofmap, bs,
1466 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1467 fn, constants,
1468 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1469 cell_info0, perms);
1470 }
1471 }
1472
1473 for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i)
1474 {
1475 auto fn = L.kernel(IntegralType::vertex, i, 0);
1476 assert(fn);
1477
1478 std::span vertices = L.domain(IntegralType::vertex, i, cell_type_idx);
1479 std::span vertices0
1480 = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx);
1481
1482 auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i});
1483
1484 assert(vertices.size() * cstride == coeffs.size());
1485
1486 impl::assemble_vertices<T>(
1487 P0, b, x_dofmap, x,
1488 md::mdspan<const std::int32_t,
1489 md::extents<std::size_t, md::dynamic_extent, 2>>(
1490 vertices.data(), vertices.size() / 2, 2),
1491 {dofs, bs,
1492 md::mdspan<const std::int32_t,
1493 md::extents<std::size_t, md::dynamic_extent, 2>>(
1494 vertices0.data(), vertices0.size() / 2, 2)},
1495 fn, constants,
1496 md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0);
1497 }
1498 }
1499}
1500
1507template <typename V, std::floating_point U,
1508 dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
1509 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
1510void assemble_vector(
1511 V&& b, const Form<T, U>& L, std::span<const T> constants,
1512 const std::map<std::pair<IntegralType, int>,
1513 std::pair<std::span<const T>, int>>& coefficients)
1514{
1515 using mdspanx3_t
1516 = md::mdspan<const scalar_value_t<T>,
1517 md::extents<std::size_t, md::dynamic_extent, 3>>;
1518
1519 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1520 assert(mesh);
1521 auto x = mesh->geometry().x();
1522 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
1523 {
1524 impl::assemble_vector(b, L, mdspanx3_t(x.data(), x.size() / 3, 3),
1525 constants, coefficients);
1526 }
1527 else
1528 {
1529 std::vector<scalar_value_t<T>> _x(x.begin(), x.end());
1530 impl::assemble_vector(b, L, mdspanx3_t(_x.data(), _x.size() / 3, 3),
1531 constants, coefficients);
1532 }
1533}
1534} // namespace dolfinx::fem::impl
Degree-of-freedom map representations and tools.
Definition DirichletBC.h:255
Functions supporting finite element method operations.
void cells(la::SparsityPattern &pattern, std::array< std::span< const std::int32_t >, 2 > cells, std::array< std::reference_wrapper< const DofMap >, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition sparsitybuild.cpp:16
Finite element method functionality.
Definition assemble_expression_impl.h:23
@ transpose
Transpose.
Definition FiniteElement.h:28
@ standard
Standard.
Definition FiniteElement.h:27
@ vertex
Vertex.
Definition Form.h:43
@ interior_facet
Interior facet.
Definition Form.h:42
@ cell
Cell.
Definition Form.h:40
@ exterior_facet
Exterior facet.
Definition Form.h:41
CellType
Cell type identifier.
Definition cell_types.h:22
int cell_num_entities(CellType type, int dim)
Definition cell_types.cpp:139