DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
assemble_vector_impl.h
1// Copyright (C) 2018-2021 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 "Constant.h"
10#include "DirichletBC.h"
11#include "DofMap.h"
12#include "Form.h"
13#include "FunctionSpace.h"
14#include "traits.h"
15#include "utils.h"
16#include <algorithm>
17#include <basix/mdspan.hpp>
18#include <cstdint>
19#include <dolfinx/common/IndexMap.h>
20#include <dolfinx/mesh/Geometry.h>
21#include <dolfinx/mesh/Mesh.h>
22#include <dolfinx/mesh/Topology.h>
23#include <functional>
24#include <memory>
25#include <optional>
26#include <span>
27#include <vector>
28
29namespace dolfinx::fem::impl
30{
31
33using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
34 const std::int32_t,
35 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
37
75template <dolfinx::scalar T, int _bs0 = -1, int _bs1 = -1>
76void _lift_bc_cells(
77 std::span<T> b, mdspan2_t x_dofmap,
78 std::span<const scalar_value_type_t<T>> x, FEkernel<T> auto kernel,
79 std::span<const std::int32_t> cells,
80 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap0,
81 fem::DofTransformKernel<T> auto P0,
82 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap1,
83 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
84 std::span<const T> coeffs, int cstride,
85 std::span<const std::uint32_t> cell_info0,
86 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
87 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha)
88{
89 if (cells.empty())
90 return;
91
92 const auto [dmap0, bs0, cells0] = dofmap0;
93 const auto [dmap1, bs1, cells1] = dofmap1;
94 assert(_bs0 < 0 or _bs0 == bs0);
95 assert(_bs1 < 0 or _bs1 == bs1);
96
97 // Data structures used in bc application
98 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * x_dofmap.extent(1));
99 std::vector<T> Ae, be;
100 assert(cells0.size() == cells.size());
101 assert(cells1.size() == cells.size());
102 for (std::size_t index = 0; index < cells.size(); ++index)
103 {
104 // Cell index in integration domain mesh, test function mesh, and trial
105 // function mesh
106 std::int32_t c = cells[index];
107 std::int32_t c0 = cells0[index];
108 std::int32_t c1 = cells1[index];
109
110 // Get dof maps for cell
111 auto dofs1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
112 dmap1, c1, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
113
114 // Check if bc is applied to cell
115 bool has_bc = false;
116 for (std::size_t j = 0; j < dofs1.size(); ++j)
117 {
118 if constexpr (_bs1 > 0)
119 {
120 for (int k = 0; k < _bs1; ++k)
121 {
122 assert(_bs1 * dofs1[j] + k < (int)bc_markers1.size());
123 if (bc_markers1[_bs1 * dofs1[j] + k])
124 {
125 has_bc = true;
126 break;
127 }
128 }
129 }
130 else
131 {
132 for (int k = 0; k < bs1; ++k)
133 {
134 assert(bs1 * dofs1[j] + k < (int)bc_markers1.size());
135 if (bc_markers1[bs1 * dofs1[j] + k])
136 {
137 has_bc = true;
138 break;
139 }
140 }
141 }
142 }
143
144 if (!has_bc)
145 continue;
146
147 // Get cell coordinates/geometry
148 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
149 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
150 for (std::size_t i = 0; i < x_dofs.size(); ++i)
151 {
152 std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
153 std::next(coordinate_dofs.begin(), 3 * i));
154 }
155
156 // Size data structure for assembly
157 auto dofs0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
158 dmap0, c0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
159
160 const int num_rows = bs0 * dofs0.size();
161 const int num_cols = bs1 * dofs1.size();
162
163 const T* coeff_array = coeffs.data() + index * cstride;
164 Ae.resize(num_rows * num_cols);
165 std::ranges::fill(Ae, 0);
166 kernel(Ae.data(), coeff_array, constants.data(), coordinate_dofs.data(),
167 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 be.resize(num_rows);
173 std::ranges::fill(be, 0);
174 for (std::size_t j = 0; j < dofs1.size(); ++j)
175 {
176 if constexpr (_bs1 > 0)
177 {
178 for (int k = 0; k < _bs1; ++k)
179 {
180 const std::int32_t jj = _bs1 * dofs1[j] + k;
181 assert(jj < (int)bc_markers1.size());
182 if (bc_markers1[jj])
183 {
184 const T bc = bc_values1[jj];
185 const T _x0 = x0.empty() ? 0 : x0[jj];
186 // const T _x0 = 0;
187 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
188 for (int m = 0; m < num_rows; ++m)
189 be[m] -= Ae[m * num_cols + _bs1 * j + k] * alpha * (bc - _x0);
190 }
191 }
192 }
193 else
194 {
195 for (int k = 0; k < bs1; ++k)
196 {
197 const std::int32_t jj = bs1 * dofs1[j] + k;
198 assert(jj < (int)bc_markers1.size());
199 if (bc_markers1[jj])
200 {
201 const T bc = bc_values1[jj];
202 const T _x0 = x0.empty() ? 0 : x0[jj];
203 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
204 for (int m = 0; m < num_rows; ++m)
205 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
206 }
207 }
208 }
209 }
210
211 for (std::size_t i = 0; i < dofs0.size(); ++i)
212 {
213 if constexpr (_bs0 > 0)
214 {
215 for (int k = 0; k < _bs0; ++k)
216 b[_bs0 * dofs0[i] + k] += be[_bs0 * i + k];
217 }
218 else
219 {
220 for (int k = 0; k < bs0; ++k)
221 b[bs0 * dofs0[i] + k] += be[bs0 * i + k];
222 }
223 }
224 }
225}
226
262template <dolfinx::scalar T, int _bs = -1>
263void _lift_bc_exterior_facets(
264 std::span<T> b, mdspan2_t x_dofmap,
265 std::span<const scalar_value_type_t<T>> x, int num_facets_per_cell,
266 FEkernel<T> auto kernel, std::span<const std::int32_t> facets,
267 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap0,
268 fem::DofTransformKernel<T> auto P0,
269 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap1,
270 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
271 std::span<const T> coeffs, int cstride,
272 std::span<const std::uint32_t> cell_info0,
273 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
274 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
275 std::span<const std::uint8_t> perms)
276{
277 if (facets.empty())
278 return;
279
280 const auto [dmap0, bs0, facets0] = dofmap0;
281 const auto [dmap1, bs1, facets1] = dofmap1;
282
283 // Data structures used in bc application
284 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * x_dofmap.extent(1));
285 std::vector<T> Ae, be;
286 assert(facets.size() % 2 == 0);
287 assert(facets0.size() == facets.size());
288 assert(facets1.size() == facets.size());
289 for (std::size_t index = 0; index < facets.size(); index += 2)
290 {
291 // Cell in integration domain mesh
292 std::int32_t cell = facets[index];
293 // Cell in test function mesh
294 std::int32_t cell0 = facets0[index];
295 // Cell in trial function mesh
296 std::int32_t cell1 = facets1[index];
297
298 // Local facet index
299 std::int32_t local_facet = facets[index + 1];
300
301 // Get dof maps for cell
302 auto dofs1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
303 dmap1, cell1, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
304
305 // Check if bc is applied to cell
306 bool has_bc = false;
307 for (std::size_t j = 0; j < dofs1.size(); ++j)
308 {
309 for (int k = 0; k < bs1; ++k)
310 {
311 if (bc_markers1[bs1 * dofs1[j] + k])
312 {
313 has_bc = true;
314 break;
315 }
316 }
317 }
318
319 if (!has_bc)
320 continue;
321
322 // Get cell coordinates/geometry
323 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
324 x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
325 for (std::size_t i = 0; i < x_dofs.size(); ++i)
326 {
327 std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
328 std::next(coordinate_dofs.begin(), 3 * i));
329 }
330
331 // Size data structure for assembly
332 auto dofs0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
333 dmap0, cell0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
334
335 const int num_rows = bs0 * dofs0.size();
336 const int num_cols = bs1 * dofs1.size();
337
338 // Permutations
339 std::uint8_t perm
340 = perms.empty() ? 0 : perms[cell * num_facets_per_cell + local_facet];
341
342 const T* coeff_array = coeffs.data() + index / 2 * cstride;
343 Ae.resize(num_rows * num_cols);
344 std::ranges::fill(Ae, 0);
345 kernel(Ae.data(), coeff_array, constants.data(), coordinate_dofs.data(),
346 &local_facet, &perm);
347 P0(Ae, cell_info0, cell0, num_cols);
348 P1T(Ae, cell_info1, cell1, num_rows);
349
350 // Size data structure for assembly
351 be.resize(num_rows);
352 std::ranges::fill(be, 0);
353 for (std::size_t j = 0; j < dofs1.size(); ++j)
354 {
355 for (int k = 0; k < bs1; ++k)
356 {
357 const std::int32_t jj = bs1 * dofs1[j] + k;
358 if (bc_markers1[jj])
359 {
360 const T bc = bc_values1[jj];
361 const T _x0 = x0.empty() ? 0 : x0[jj];
362 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
363 for (int m = 0; m < num_rows; ++m)
364 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
365 }
366 }
367 }
368
369 for (std::size_t i = 0; i < dofs0.size(); ++i)
370 for (int k = 0; k < bs0; ++k)
371 b[bs0 * dofs0[i] + k] += be[bs0 * i + k];
372 }
373}
374
411template <dolfinx::scalar T, int _bs = -1>
412void _lift_bc_interior_facets(
413 std::span<T> b, mdspan2_t x_dofmap,
414 std::span<const scalar_value_type_t<T>> x, int num_facets_per_cell,
415 FEkernel<T> auto kernel, std::span<const std::int32_t> facets,
416 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap0,
417 fem::DofTransformKernel<T> auto P0,
418 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap1,
419 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
420 std::span<const T> coeffs, int cstride,
421 std::span<const std::uint32_t> cell_info0,
422 std::span<const std::uint32_t> cell_info1,
423 std::span<const std::uint8_t> perms, std::span<const T> bc_values1,
424 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha)
425{
426 if (facets.empty())
427 return;
428
429 const auto [dmap0, bs0, facets0] = dofmap0;
430 const auto [dmap1, bs1, facets1] = dofmap1;
431
432 // Data structures used in assembly
433 using X = scalar_value_type_t<T>;
434 std::vector<X> coordinate_dofs(2 * x_dofmap.extent(1) * 3);
435 std::span<X> cdofs0(coordinate_dofs.data(), x_dofmap.extent(1) * 3);
436 std::span<X> cdofs1(coordinate_dofs.data() + x_dofmap.extent(1) * 3,
437 x_dofmap.extent(1) * 3);
438 std::vector<T> Ae, be;
439
440 // Temporaries for joint dofmaps
441 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
442 assert(facets.size() % 4 == 0);
443
444 const int num_dofs0 = dmap0.extent(1);
445 const int num_dofs1 = dmap1.extent(1);
446 assert(facets0.size() == facets.size());
447 assert(facets1.size() == facets.size());
448 for (std::size_t index = 0; index < facets.size(); index += 4)
449 {
450 // Cells in integration domain, test function domain and trial
451 // function domain meshes
452 std::array cells{facets[index], facets[index + 2]};
453 std::array cells0{facets0[index], facets0[index + 2]};
454 std::array cells1{facets1[index], facets1[index + 2]};
455
456 // Local facet indices
457 std::array<std::int32_t, 2> local_facet
458 = {facets[index + 1], facets[index + 3]};
459
460 // Get cell geometry
461 auto x_dofs0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
462 x_dofmap, cells[0], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
463 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
464 {
465 std::copy_n(std::next(x.begin(), 3 * x_dofs0[i]), 3,
466 std::next(cdofs0.begin(), 3 * i));
467 }
468 auto x_dofs1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
469 x_dofmap, cells[1], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
470 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
471 {
472 std::copy_n(std::next(x.begin(), 3 * x_dofs1[i]), 3,
473 std::next(cdofs1.begin(), 3 * i));
474 }
475
476 // Get dof maps for cells and pack
477 auto dmap0_cell0
478 = std::span(dmap0.data_handle() + cells0[0] * num_dofs0, num_dofs0);
479 auto dmap0_cell1
480 = std::span(dmap0.data_handle() + cells0[1] * num_dofs0, num_dofs0);
481
482 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
483 std::ranges::copy(dmap0_cell0, dmapjoint0.begin());
484 std::ranges::copy(dmap0_cell1,
485 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
486
487 auto dmap1_cell0
488 = std::span(dmap1.data_handle() + cells1[0] * num_dofs1, num_dofs1);
489 auto dmap1_cell1
490 = std::span(dmap1.data_handle() + cells1[1] * num_dofs1, num_dofs1);
491
492 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
493 std::ranges::copy(dmap1_cell0, dmapjoint1.begin());
494 std::ranges::copy(dmap1_cell1,
495 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
496
497 // Check if bc is applied to cell0
498 bool has_bc = false;
499 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
500 {
501 for (int k = 0; k < bs1; ++k)
502 {
503 if (bc_markers1[bs1 * dmap1_cell0[j] + k])
504 {
505 has_bc = true;
506 break;
507 }
508 }
509 }
510
511 // Check if bc is applied to cell1
512 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
513 {
514 for (int k = 0; k < bs1; ++k)
515 {
516 if (bc_markers1[bs1 * dmap1_cell1[j] + k])
517 {
518 has_bc = true;
519 break;
520 }
521 }
522 }
523
524 if (!has_bc)
525 continue;
526
527 const int num_rows = bs0 * dmapjoint0.size();
528 const int num_cols = bs1 * dmapjoint1.size();
529
530 // Tabulate tensor
531 Ae.resize(num_rows * num_cols);
532 std::ranges::fill(Ae, 0);
533 std::array perm
534 = perms.empty()
535 ? std::array<std::uint8_t, 2>{0, 0}
536 : std::array{
537 perms[cells[0] * num_facets_per_cell + local_facet[0]],
538 perms[cells[1] * num_facets_per_cell + local_facet[1]]};
539 kernel(Ae.data(), coeffs.data() + index / 2 * cstride, constants.data(),
540 coordinate_dofs.data(), local_facet.data(), perm.data());
541
542 std::span<T> _Ae(Ae);
543 std::span<T> sub_Ae0 = _Ae.subspan(bs0 * dmap0_cell0.size() * num_cols,
544 bs0 * dmap0_cell1.size() * num_cols);
545
546 P0(_Ae, cell_info0, cells0[0], num_cols);
547 P0(sub_Ae0, cell_info0, cells0[1], num_cols);
548 P1T(_Ae, cell_info1, cells1[0], num_rows);
549
550 for (int row = 0; row < num_rows; ++row)
551 {
552 // DOFs for dmap1 and cell1 are not stored contiguously in
553 // the block matrix, so each row needs a separate span access
554 std::span<T> sub_Ae1 = _Ae.subspan(
555 row * num_cols + bs1 * dmap1_cell0.size(), bs1 * dmap1_cell1.size());
556 P1T(sub_Ae1, cell_info1, cells1[1], 1);
557 }
558
559 be.resize(num_rows);
560 std::ranges::fill(be, 0);
561
562 // Compute b = b - A*b for cell0
563 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
564 {
565 for (int k = 0; k < bs1; ++k)
566 {
567 const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
568 if (bc_markers1[jj])
569 {
570 const T bc = bc_values1[jj];
571 const T _x0 = x0.empty() ? 0 : x0[jj];
572 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
573 for (int m = 0; m < num_rows; ++m)
574 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
575 }
576 }
577 }
578
579 // Compute b = b - A*b for cell1
580 const int offset = bs1 * dmap1_cell0.size();
581 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
582 {
583 for (int k = 0; k < bs1; ++k)
584 {
585 const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
586 if (bc_markers1[jj])
587 {
588 const T bc = bc_values1[jj];
589 const T _x0 = x0.empty() ? 0 : x0[jj];
590 // be -= Ae.col(offset + bs1 * j + k) * alpha * (bc - x0[jj]);
591 for (int m = 0; m < num_rows; ++m)
592 {
593 be[m]
594 -= Ae[m * num_cols + offset + bs1 * j + k] * alpha * (bc - _x0);
595 }
596 }
597 }
598 }
599
600 for (std::size_t i = 0; i < dmap0_cell0.size(); ++i)
601 for (int k = 0; k < bs0; ++k)
602 b[bs0 * dmap0_cell0[i] + k] += be[bs0 * i + k];
603
604 const int offset_be = bs0 * dmap0_cell0.size();
605 for (std::size_t i = 0; i < dmap0_cell1.size(); ++i)
606 for (int k = 0; k < bs0; ++k)
607 b[bs0 * dmap0_cell1[i] + k] += be[offset_be + bs0 * i + k];
608 }
609}
610
633template <dolfinx::scalar T, int _bs = -1>
634void assemble_cells(
635 fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
636 std::span<const scalar_value_type_t<T>> x,
637 std::span<const std::int32_t> cells,
638 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap,
639 FEkernel<T> auto kernel, std::span<const T> constants,
640 std::span<const T> coeffs, int cstride,
641 std::span<const std::uint32_t> cell_info0)
642{
643 if (cells.empty())
644 return;
645
646 const auto [dmap, bs, cells0] = dofmap;
647 assert(_bs < 0 or _bs == bs);
648
649 // Create data structures used in assembly
650 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * x_dofmap.extent(1));
651 std::vector<T> be(bs * dmap.extent(1));
652 std::span<T> _be(be);
653
654 // Iterate over active cells
655 for (std::size_t index = 0; index < cells.size(); ++index)
656 {
657 // Integration domain celland test function cell
658 std::int32_t c = cells[index];
659 std::int32_t c0 = cells0[index];
660
661 // Get cell coordinates/geometry
662 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
663 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
664 for (std::size_t i = 0; i < x_dofs.size(); ++i)
665 {
666 std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
667 std::next(coordinate_dofs.begin(), 3 * i));
668 }
669
670 // Tabulate vector for cell
671 std::ranges::fill(be, 0);
672 kernel(be.data(), coeffs.data() + index * cstride, constants.data(),
673 coordinate_dofs.data(), nullptr, nullptr);
674 P0(_be, cell_info0, c0, 1);
675
676 // Scatter cell vector to 'global' vector array
677 auto dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
678 dmap, c0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
679 if constexpr (_bs > 0)
680 {
681 for (std::size_t i = 0; i < dofs.size(); ++i)
682 for (int k = 0; k < _bs; ++k)
683 b[_bs * dofs[i] + k] += be[_bs * i + k];
684 }
685 else
686 {
687 for (std::size_t i = 0; i < dofs.size(); ++i)
688 for (int k = 0; k < bs; ++k)
689 b[bs * dofs[i] + k] += be[bs * i + k];
690 }
691 }
692}
693
719template <dolfinx::scalar T, int _bs = -1>
720void assemble_exterior_facets(
721 fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
722 std::span<const scalar_value_type_t<T>> x, int num_facets_per_cell,
723 std::span<const std::int32_t> facets,
724 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap,
725 FEkernel<T> auto fn, std::span<const T> constants,
726 std::span<const T> coeffs, int cstride,
727 std::span<const std::uint32_t> cell_info0,
728 std::span<const std::uint8_t> perms)
729{
730 if (facets.empty())
731 return;
732
733 const auto [dmap, bs, facets0] = dofmap;
734 assert(_bs < 0 or _bs == bs);
735
736 // FIXME: Add proper interface for num_dofs
737 // Create data structures used in assembly
738 const int num_dofs = dmap.extent(1);
739 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * x_dofmap.extent(1));
740 std::vector<T> be(bs * num_dofs);
741 std::span<T> _be(be);
742 assert(facets.size() % 2 == 0);
743 assert(facets0.size() == facets.size());
744 for (std::size_t index = 0; index < facets.size(); index += 2)
745 {
746 // Cell in the integration domain, local facet index relative to the
747 // integration domain cell, and cell in the test function mesh
748 std::int32_t cell = facets[index];
749 std::int32_t local_facet = facets[index + 1];
750 std::int32_t cell0 = facets0[index];
751
752 // Get cell coordinates/geometry
753 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
754 x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
755 for (std::size_t i = 0; i < x_dofs.size(); ++i)
756 {
757 std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
758 std::next(coordinate_dofs.begin(), 3 * i));
759 }
760
761 // Permutations
762 std::uint8_t perm
763 = perms.empty() ? 0 : perms[cell * num_facets_per_cell + local_facet];
764
765 // Tabulate element vector
766 std::ranges::fill(be, 0);
767 fn(be.data(), coeffs.data() + index / 2 * cstride, constants.data(),
768 coordinate_dofs.data(), &local_facet, &perm);
769
770 P0(_be, cell_info0, cell0, 1);
771
772 // Add element vector to global vector
773 auto dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
774 dmap, cell0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
775 if constexpr (_bs > 0)
776 {
777 for (std::size_t i = 0; i < dofs.size(); ++i)
778 for (int k = 0; k < _bs; ++k)
779 b[_bs * dofs[i] + k] += be[_bs * i + k];
780 }
781 else
782 {
783 for (std::size_t i = 0; i < dofs.size(); ++i)
784 for (int k = 0; k < bs; ++k)
785 b[bs * dofs[i] + k] += be[bs * i + k];
786 }
787 }
788}
789
815template <dolfinx::scalar T, int _bs = -1>
816void assemble_interior_facets(
817 fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
818 std::span<const scalar_value_type_t<T>> x, int num_facets_per_cell,
819 std::span<const std::int32_t> facets,
820 std::tuple<const DofMap&, int, std::span<const std::int32_t>> dofmap,
821 FEkernel<T> auto fn, std::span<const T> constants,
822 std::span<const T> coeffs, int cstride,
823 std::span<const std::uint32_t> cell_info0,
824 std::span<const std::uint8_t> perms)
825{
826 if (facets.empty())
827 return;
828
829 const auto [dmap, bs, facets0] = dofmap;
830 assert(_bs < 0 or _bs == bs);
831
832 // Create data structures used in assembly
833 using X = scalar_value_type_t<T>;
834 std::vector<X> coordinate_dofs(2 * x_dofmap.extent(1) * 3);
835 std::span<X> cdofs0(coordinate_dofs.data(), x_dofmap.extent(1) * 3);
836 std::span<X> cdofs1(coordinate_dofs.data() + x_dofmap.extent(1) * 3,
837 x_dofmap.extent(1) * 3);
838 std::vector<T> be;
839
840 assert(facets.size() % 4 == 0);
841 assert(facets0.size() == facets.size());
842 for (std::size_t index = 0; index < facets.size(); index += 4)
843 {
844 // Cells in integration domain and test function domain meshes
845 std::array<std::int32_t, 2> cells{facets[index], facets[index + 2]};
846 std::array<std::int32_t, 2> cells0{facets0[index], facets0[index + 2]};
847
848 // Local facet indices
849 std::array<std::int32_t, 2> local_facet{facets[index + 1],
850 facets[index + 3]};
851
852 // Get cell geometry
853 auto x_dofs0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
854 x_dofmap, cells[0], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
855 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
856 {
857 std::copy_n(std::next(x.begin(), 3 * x_dofs0[i]), 3,
858 std::next(cdofs0.begin(), 3 * i));
859 }
860 auto x_dofs1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
861 x_dofmap, cells[1], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
862 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
863 {
864 std::copy_n(std::next(x.begin(), 3 * x_dofs1[i]), 3,
865 std::next(cdofs1.begin(), 3 * i));
866 }
867
868 // Get dofmaps for cells
869 std::span<const std::int32_t> dmap0 = dmap.cell_dofs(cells0[0]);
870 std::span<const std::int32_t> dmap1 = dmap.cell_dofs(cells0[1]);
871
872 // Tabulate element vector
873 be.resize(bs * (dmap0.size() + dmap1.size()));
874 std::ranges::fill(be, 0);
875 std::array perm
876 = perms.empty()
877 ? std::array<std::uint8_t, 2>{0, 0}
878 : std::array{
879 perms[cells[0] * num_facets_per_cell + local_facet[0]],
880 perms[cells[1] * num_facets_per_cell + local_facet[1]]};
881 fn(be.data(), coeffs.data() + index / 2 * cstride, constants.data(),
882 coordinate_dofs.data(), local_facet.data(), perm.data());
883
884 std::span<T> _be(be);
885 std::span<T> sub_be = _be.subspan(bs * dmap0.size(), bs * dmap1.size());
886
887 P0(be, cell_info0, cells0[0], 1);
888 P0(sub_be, cell_info0, cells0[1], 1);
889
890 // Add element vector to global vector
891 if constexpr (_bs > 0)
892 {
893 for (std::size_t i = 0; i < dmap0.size(); ++i)
894 for (int k = 0; k < _bs; ++k)
895 b[_bs * dmap0[i] + k] += be[_bs * i + k];
896 for (std::size_t i = 0; i < dmap1.size(); ++i)
897 for (int k = 0; k < _bs; ++k)
898 b[_bs * dmap1[i] + k] += be[_bs * (i + dmap0.size()) + k];
899 }
900 else
901 {
902 for (std::size_t i = 0; i < dmap0.size(); ++i)
903 for (int k = 0; k < bs; ++k)
904 b[bs * dmap0[i] + k] += be[bs * i + k];
905 for (std::size_t i = 0; i < dmap1.size(); ++i)
906 for (int k = 0; k < bs; ++k)
907 b[bs * dmap1[i] + k] += be[bs * (i + dmap0.size()) + k];
908 }
909 }
910}
911
928template <dolfinx::scalar T, std::floating_point U>
929void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
930 std::span<const scalar_value_type_t<T>> x,
931 std::span<const T> constants,
932 const std::map<std::pair<IntegralType, int>,
933 std::pair<std::span<const T>, int>>& coefficients,
934 std::span<const T> bc_values1,
935 std::span<const std::int8_t> bc_markers1, std::span<const T> x0,
936 T alpha)
937{
938 // Integration domain mesh
939 std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
940 assert(mesh);
941 // Test function mesh
942 auto mesh0 = a.function_spaces().at(0)->mesh();
943 assert(mesh0);
944 // Trial function mesh
945 auto mesh1 = a.function_spaces().at(1)->mesh();
946 assert(mesh1);
947
948 // Get dofmap for columns and rows of a
949 assert(a.function_spaces().at(0));
950 assert(a.function_spaces().at(1));
951 auto dofmap0 = a.function_spaces()[0]->dofmap()->map();
952 const int bs0 = a.function_spaces()[0]->dofmap()->bs();
953 auto element0 = a.function_spaces()[0]->element();
954 auto dofmap1 = a.function_spaces()[1]->dofmap()->map();
955 const int bs1 = a.function_spaces()[1]->dofmap()->bs();
956 auto element1 = a.function_spaces()[1]->element();
957 assert(element0);
958
959 std::span<const std::uint32_t> cell_info0;
960 std::span<const std::uint32_t> cell_info1;
961 // TODO: Check for each element instead
962 if (element0->needs_dof_transformations()
963 or element1->needs_dof_transformations() or a.needs_facet_permutations())
964 {
965 mesh0->topology_mutable()->create_entity_permutations();
966 mesh1->topology_mutable()->create_entity_permutations();
967 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
968 cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
969 }
970
971 fem::DofTransformKernel<T> auto P0
972 = element0->template dof_transformation_fn<T>(doftransform::standard);
973 fem::DofTransformKernel<T> auto P1T
974 = element1->template dof_transformation_right_fn<T>(
976
977 for (int i : a.integral_ids(IntegralType::cell))
978 {
979 auto kernel = a.kernel(IntegralType::cell, i);
980 assert(kernel);
981 auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
982 std::span<const std::int32_t> cells = a.domain(IntegralType::cell, i);
983 if (bs0 == 1 and bs1 == 1)
984 {
985 _lift_bc_cells<T, 1, 1>(
986 b, x_dofmap, x, kernel, cells,
987 {dofmap0, bs0, a.domain(IntegralType::cell, i, *mesh0)}, P0,
988 {dofmap1, bs1, a.domain(IntegralType::cell, i, *mesh1)}, P1T,
989 constants, coeffs, cstride, cell_info0, cell_info1, bc_values1,
990 bc_markers1, x0, alpha);
991 }
992 else if (bs0 == 3 and bs1 == 3)
993 {
994 _lift_bc_cells<T, 3, 3>(
995 b, x_dofmap, x, kernel, cells,
996 {dofmap0, bs0, a.domain(IntegralType::cell, i, *mesh0)}, P0,
997 {dofmap1, bs1, a.domain(IntegralType::cell, i, *mesh1)}, P1T,
998 constants, coeffs, cstride, cell_info0, cell_info1, bc_values1,
999 bc_markers1, x0, alpha);
1000 }
1001 else
1002 {
1003 _lift_bc_cells(b, x_dofmap, x, kernel, cells,
1004 {dofmap0, bs0, a.domain(IntegralType::cell, i, *mesh0)},
1005 P0,
1006 {dofmap1, bs1, a.domain(IntegralType::cell, i, *mesh1)},
1007 P1T, constants, coeffs, cstride, cell_info0, cell_info1,
1008 bc_values1, bc_markers1, x0, alpha);
1009 }
1010 }
1011
1012 std::span<const std::uint8_t> perms;
1013 if (a.needs_facet_permutations())
1014 {
1015 mesh->topology_mutable()->create_entity_permutations();
1016 perms = std::span(mesh->topology()->get_facet_permutations());
1017 }
1018
1019 mesh::CellType cell_type = mesh->topology()->cell_type();
1020 int num_facets_per_cell
1021 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
1022 for (int i : a.integral_ids(IntegralType::exterior_facet))
1023 {
1024 auto kernel = a.kernel(IntegralType::exterior_facet, i);
1025 assert(kernel);
1026 auto& [coeffs, cstride]
1027 = coefficients.at({IntegralType::exterior_facet, i});
1028 _lift_bc_exterior_facets(
1029 b, x_dofmap, x, num_facets_per_cell, kernel,
1030 a.domain(IntegralType::exterior_facet, i),
1031 {dofmap0, bs0, a.domain(IntegralType::exterior_facet, i, *mesh0)}, P0,
1032 {dofmap1, bs1, a.domain(IntegralType::exterior_facet, i, *mesh1)}, P1T,
1033 constants, coeffs, cstride, cell_info0, cell_info1, bc_values1,
1034 bc_markers1, x0, alpha, perms);
1035 }
1036
1037 for (int i : a.integral_ids(IntegralType::interior_facet))
1038 {
1039 auto kernel = a.kernel(IntegralType::interior_facet, i);
1040 assert(kernel);
1041 auto& [coeffs, cstride]
1042 = coefficients.at({IntegralType::interior_facet, i});
1043 _lift_bc_interior_facets(
1044 b, x_dofmap, x, num_facets_per_cell, kernel,
1045 a.domain(IntegralType::interior_facet, i),
1046 {dofmap0, bs0, a.domain(IntegralType::interior_facet, i, *mesh0)}, P0,
1047 {dofmap1, bs1, a.domain(IntegralType::interior_facet, i, *mesh1)}, P1T,
1048 constants, coeffs, cstride, cell_info0, cell_info1, perms, bc_values1,
1049 bc_markers1, x0, alpha);
1050 }
1051}
1052
1073template <dolfinx::scalar T, std::floating_point U>
1074void apply_lifting(
1075 std::span<T> b, const std::vector<std::shared_ptr<const Form<T, U>>> a,
1076 const std::vector<std::span<const T>>& constants,
1077 const std::vector<std::map<std::pair<IntegralType, int>,
1078 std::pair<std::span<const T>, int>>>& coeffs,
1079 const std::vector<std::vector<std::shared_ptr<const DirichletBC<T, U>>>>&
1080 bcs1,
1081 const std::vector<std::span<const T>>& x0, T alpha)
1082{
1083 if (!x0.empty() and x0.size() != a.size())
1084 {
1085 throw std::runtime_error(
1086 "Mismatch in size between x0 and bilinear form in assembler.");
1087 }
1088
1089 if (a.size() != bcs1.size())
1090 {
1091 throw std::runtime_error(
1092 "Mismatch in size between a and bcs in assembler.");
1093 }
1094
1095 for (std::size_t j = 0; j < a.size(); ++j)
1096 {
1097 std::vector<std::int8_t> bc_markers1;
1098 std::vector<T> bc_values1;
1099 if (a[j] and !bcs1[j].empty())
1100 {
1101 // Extract data from mesh
1102 std::shared_ptr<const mesh::Mesh<U>> mesh = a[j]->mesh();
1103 if (!mesh)
1104 throw std::runtime_error("Unable to extract a mesh.");
1105 mdspan2_t x_dofmap = mesh->geometry().dofmap();
1106 auto x = mesh->geometry().x();
1107
1108 assert(a[j]->function_spaces().at(0));
1109
1110 auto V1 = a[j]->function_spaces()[1];
1111 assert(V1);
1112 auto map1 = V1->dofmap()->index_map;
1113 const int bs1 = V1->dofmap()->index_map_bs();
1114 assert(map1);
1115 const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
1116 bc_markers1.assign(crange, false);
1117 bc_values1.assign(crange, 0);
1118 for (const std::shared_ptr<const DirichletBC<T, U>>& bc : bcs1[j])
1119 {
1120 bc->mark_dofs(bc_markers1);
1121 bc->set(bc_values1, std::nullopt, 1);
1122 }
1123
1124 if (!x0.empty())
1125 {
1126 lift_bc<T>(b, *a[j], x_dofmap, x, constants[j], coeffs[j], bc_values1,
1127 bc_markers1, x0[j], alpha);
1128 }
1129 else
1130 {
1131 lift_bc<T>(b, *a[j], x_dofmap, x, constants[j], coeffs[j], bc_values1,
1132 bc_markers1, std::span<const T>(), alpha);
1133 }
1134 }
1135 }
1136}
1137
1146template <dolfinx::scalar T, std::floating_point U>
1147void assemble_vector(
1148 std::span<T> b, const Form<T, U>& L, mdspan2_t x_dofmap,
1149 std::span<const scalar_value_type_t<T>> x, std::span<const T> constants,
1150 const std::map<std::pair<IntegralType, int>,
1151 std::pair<std::span<const T>, int>>& coefficients)
1152{
1153 // Integration domain mesh
1154 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1155 assert(mesh);
1156
1157 // Test function mesh
1158 auto mesh0 = L.function_spaces().at(0)->mesh();
1159 assert(mesh0);
1160
1161 // Get dofmap data
1162 assert(L.function_spaces().at(0));
1163 auto element = L.function_spaces().at(0)->element();
1164 assert(element);
1165 std::shared_ptr<const fem::DofMap> dofmap
1166 = L.function_spaces().at(0)->dofmap();
1167 assert(dofmap);
1168 auto dofs = dofmap->map();
1169 const int bs = dofmap->bs();
1170
1171 fem::DofTransformKernel<T> auto P0
1172 = element->template dof_transformation_fn<T>(doftransform::standard);
1173
1174 std::span<const std::uint32_t> cell_info0;
1175 if (element->needs_dof_transformations() or L.needs_facet_permutations())
1176 {
1177 mesh0->topology_mutable()->create_entity_permutations();
1178 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
1179 }
1180
1181 for (int i : L.integral_ids(IntegralType::cell))
1182 {
1183 auto fn = L.kernel(IntegralType::cell, i);
1184 assert(fn);
1185 auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
1186 std::span<const std::int32_t> cells = L.domain(IntegralType::cell, i);
1187 if (bs == 1)
1188 {
1189 impl::assemble_cells<T, 1>(
1190 P0, b, x_dofmap, x, cells,
1191 {dofs, bs, L.domain(IntegralType::cell, i, *mesh0)}, fn, constants,
1192 coeffs, cstride, cell_info0);
1193 }
1194 else if (bs == 3)
1195 {
1196 impl::assemble_cells<T, 3>(
1197 P0, b, x_dofmap, x, cells,
1198 {dofs, bs, L.domain(IntegralType::cell, i, *mesh0)}, fn, constants,
1199 coeffs, cstride, cell_info0);
1200 }
1201 else
1202 {
1203 impl::assemble_cells(P0, b, x_dofmap, x, cells,
1204 {dofs, bs, L.domain(IntegralType::cell, i, *mesh0)},
1205 fn, constants, coeffs, cstride, cell_info0);
1206 }
1207 }
1208
1209 std::span<const std::uint8_t> perms;
1210 if (L.needs_facet_permutations())
1211 {
1212 mesh->topology_mutable()->create_entity_permutations();
1213 perms = std::span(mesh->topology()->get_facet_permutations());
1214 }
1215
1216 mesh::CellType cell_type = mesh->topology()->cell_type();
1217 int num_facets_per_cell
1218 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
1219 for (int i : L.integral_ids(IntegralType::exterior_facet))
1220 {
1221 auto fn = L.kernel(IntegralType::exterior_facet, i);
1222 assert(fn);
1223 auto& [coeffs, cstride]
1224 = coefficients.at({IntegralType::exterior_facet, i});
1225 std::span<const std::int32_t> facets
1226 = L.domain(IntegralType::exterior_facet, i);
1227 if (bs == 1)
1228 {
1229 impl::assemble_exterior_facets<T, 1>(
1230 P0, b, x_dofmap, x, num_facets_per_cell, facets,
1231 {dofs, bs, L.domain(IntegralType::exterior_facet, i, *mesh0)}, fn,
1232 constants, coeffs, cstride, cell_info0, perms);
1233 }
1234 else if (bs == 3)
1235 {
1236 impl::assemble_exterior_facets<T, 3>(
1237 P0, b, x_dofmap, x, num_facets_per_cell, facets,
1238 {dofs, bs, L.domain(IntegralType::exterior_facet, i, *mesh0)}, fn,
1239 constants, coeffs, cstride, cell_info0, perms);
1240 }
1241 else
1242 {
1243 impl::assemble_exterior_facets(
1244 P0, b, x_dofmap, x, num_facets_per_cell, facets,
1245 {dofs, bs, L.domain(IntegralType::exterior_facet, i, *mesh0)}, fn,
1246 constants, coeffs, cstride, cell_info0, perms);
1247 }
1248 }
1249
1250 for (int i : L.integral_ids(IntegralType::interior_facet))
1251 {
1252 auto fn = L.kernel(IntegralType::interior_facet, i);
1253 assert(fn);
1254 auto& [coeffs, cstride]
1255 = coefficients.at({IntegralType::interior_facet, i});
1256 std::span<const std::int32_t> facets
1257 = L.domain(IntegralType::interior_facet, i);
1258 if (bs == 1)
1259 {
1260 impl::assemble_interior_facets<T, 1>(
1261 P0, b, x_dofmap, x, num_facets_per_cell, facets,
1262 {*dofmap, bs, L.domain(IntegralType::interior_facet, i, *mesh0)}, fn,
1263 constants, coeffs, cstride, cell_info0, perms);
1264 }
1265 else if (bs == 3)
1266 {
1267 impl::assemble_interior_facets<T, 3>(
1268 P0, b, x_dofmap, x, num_facets_per_cell, facets,
1269 {*dofmap, bs, L.domain(IntegralType::interior_facet, i, *mesh0)}, fn,
1270 constants, coeffs, cstride, cell_info0, perms);
1271 }
1272 else
1273 {
1274 impl::assemble_interior_facets(
1275 P0, b, x_dofmap, x, num_facets_per_cell, facets,
1276 {*dofmap, bs, L.domain(IntegralType::interior_facet, i, *mesh0)}, fn,
1277 constants, coeffs, cstride, cell_info0, perms);
1278 }
1279 }
1280}
1281
1288template <dolfinx::scalar T, std::floating_point U>
1289void assemble_vector(
1290 std::span<T> b, const Form<T, U>& L, std::span<const T> constants,
1291 const std::map<std::pair<IntegralType, int>,
1292 std::pair<std::span<const T>, int>>& coefficients)
1293{
1294 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1295 assert(mesh);
1296 if constexpr (std::is_same_v<U, scalar_value_type_t<T>>)
1297 assemble_vector(b, L, mesh->geometry().dofmap(), mesh->geometry().x(),
1298 constants, coefficients);
1299 else
1300 {
1301 auto x = mesh->geometry().x();
1302 std::vector<scalar_value_type_t<T>> _x(x.begin(), x.end());
1303 assemble_vector(b, L, mesh->geometry().dofmap(), _x, constants,
1304 coefficients);
1305 }
1306}
1307} // namespace dolfinx::fem::impl
Degree-of-freedom map representations and tools.
Definition types.h:20
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
@ interior_facet
Interior facet.
@ exterior_facet
Exterior facet.
int cell_num_entities(CellType type, int dim)
Definition cell_types.cpp:139
CellType
Cell type identifier.
Definition cell_types.h:22