DOLFINx 0.10.0.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 "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 <dolfinx::scalar T, int _bs0 = -1, int _bs1 = -1>
80void _lift_bc_cells(
81 std::span<T> b, mdspan2_t x_dofmap,
82 md::mdspan<const scalar_value_t<T>,
83 md::extents<std::size_t, md::dynamic_extent, 3>>
84 x,
85 FEkernel<T> auto kernel, std::span<const std::int32_t> cells,
86 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap0,
87 fem::DofTransformKernel<T> auto P0,
88 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap1,
89 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
90 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
91 std::span<const std::uint32_t> cell_info0,
92 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
93 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha)
94{
95 if (cells.empty())
96 return;
97
98 const auto [dmap0, bs0, cells0] = dofmap0;
99 const auto [dmap1, bs1, cells1] = dofmap1;
100 assert(_bs0 < 0 or _bs0 == bs0);
101 assert(_bs1 < 0 or _bs1 == bs1);
102
103 // Data structures used in bc application
104 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
105 std::vector<T> Ae, be;
106 assert(cells0.size() == cells.size());
107 assert(cells1.size() == cells.size());
108 for (std::size_t index = 0; index < cells.size(); ++index)
109 {
110 // Cell index in integration domain mesh, test function mesh, and trial
111 // function mesh
112 std::int32_t c = cells[index];
113 std::int32_t c0 = cells0[index];
114 std::int32_t c1 = cells1[index];
115
116 // Get dof maps for cell
117 auto dofs1 = md::submdspan(dmap1, c1, md::full_extent);
118
119 // Check if bc is applied to cell
120 bool has_bc = false;
121 for (std::size_t j = 0; j < dofs1.size(); ++j)
122 {
123 if constexpr (_bs1 > 0)
124 {
125 for (int k = 0; k < _bs1; ++k)
126 {
127 assert(_bs1 * dofs1[j] + k < (int)bc_markers1.size());
128 if (bc_markers1[_bs1 * dofs1[j] + k])
129 {
130 has_bc = true;
131 break;
132 }
133 }
134 }
135 else
136 {
137 for (int k = 0; k < bs1; ++k)
138 {
139 assert(bs1 * dofs1[j] + k < (int)bc_markers1.size());
140 if (bc_markers1[bs1 * dofs1[j] + k])
141 {
142 has_bc = true;
143 break;
144 }
145 }
146 }
147 }
148
149 if (!has_bc)
150 continue;
151
152 // Get cell coordinates/geometry
153 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
154 for (std::size_t i = 0; i < x_dofs.size(); ++i)
155 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
156
157 // Size data structure for assembly
158 auto dofs0 = md::submdspan(dmap0, c0, md::full_extent);
159
160 const int num_rows = bs0 * dofs0.size();
161 const int num_cols = bs1 * dofs1.size();
162
163 Ae.resize(num_rows * num_cols);
164 std::ranges::fill(Ae, 0);
165 kernel(Ae.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
166 nullptr, nullptr, nullptr);
167 P0(Ae, cell_info0, c0, num_cols);
168 P1T(Ae, cell_info1, c1, num_rows);
169
170 // Size data structure for assembly
171 be.resize(num_rows);
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 <dolfinx::scalar T>
264void _lift_bc_exterior_facets(
265 std::span<T> b, mdspan2_t x_dofmap,
266 md::mdspan<const scalar_value_t<T>,
267 md::extents<std::size_t, md::dynamic_extent, 3>>
268 x,
269 FEkernel<T> auto kernel,
270 md::mdspan<const std::int32_t,
271 md::extents<std::size_t, md::dynamic_extent, 2>>
272 facets,
273 std::tuple<mdspan2_t, int,
274 md::mdspan<const std::int32_t,
275 md::extents<std::size_t, md::dynamic_extent, 2>>>
276 dofmap0,
277 fem::DofTransformKernel<T> auto P0,
278 std::tuple<mdspan2_t, int,
279 md::mdspan<const std::int32_t,
280 md::extents<std::size_t, md::dynamic_extent, 2>>>
281 dofmap1,
282 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
283 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
284 std::span<const std::uint32_t> cell_info0,
285 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
286 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
287 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
288{
289 if (facets.empty())
290 return;
291
292 const auto [dmap0, bs0, facets0] = dofmap0;
293 const auto [dmap1, bs1, facets1] = dofmap1;
294
295 // Data structures used in bc application
296 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
297 std::vector<T> Ae, be;
298 assert(facets0.size() == facets.size());
299 assert(facets1.size() == facets.size());
300 for (std::size_t index = 0; index < facets.extent(0); ++index)
301 {
302 // Cell in integration domain, test function and trial function
303 // meshes
304 std::int32_t cell = facets(index, 0);
305 std::int32_t cell0 = facets0(index, 0);
306 std::int32_t cell1 = facets1(index, 0);
307
308 // Local facet index
309 std::int32_t local_facet = facets(index, 1);
310
311 // Get dof maps for cell
312 auto dofs1 = md::submdspan(dmap1, cell1, md::full_extent);
313
314 // Check if bc is applied to cell
315 bool has_bc = false;
316 for (std::size_t j = 0; j < dofs1.size(); ++j)
317 {
318 for (int k = 0; k < bs1; ++k)
319 {
320 if (bc_markers1[bs1 * dofs1[j] + k])
321 {
322 has_bc = true;
323 break;
324 }
325 }
326 }
327
328 if (!has_bc)
329 continue;
330
331 // Get cell coordinates/geometry
332 auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
333 for (std::size_t i = 0; i < x_dofs.size(); ++i)
334 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
335
336 // Size data structure for assembly
337 auto dofs0 = md::submdspan(dmap0, cell0, md::full_extent);
338
339 const int num_rows = bs0 * dofs0.size();
340 const int num_cols = bs1 * dofs1.size();
341
342 // Permutations
343 std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
344
345 Ae.resize(num_rows * num_cols);
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 be.resize(num_rows);
354 std::ranges::fill(be, 0);
355 for (std::size_t j = 0; j < dofs1.size(); ++j)
356 {
357 for (int k = 0; k < bs1; ++k)
358 {
359 const std::int32_t jj = bs1 * dofs1[j] + k;
360 if (bc_markers1[jj])
361 {
362 const T bc = bc_values1[jj];
363 const T _x0 = x0.empty() ? 0 : x0[jj];
364 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
365 for (int m = 0; m < num_rows; ++m)
366 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
367 }
368 }
369 }
370
371 for (std::size_t i = 0; i < dofs0.size(); ++i)
372 for (int k = 0; k < bs0; ++k)
373 b[bs0 * dofs0[i] + k] += be[bs0 * i + k];
374 }
375}
376
417template <dolfinx::scalar T>
418void _lift_bc_interior_facets(
419 std::span<T> b, mdspan2_t x_dofmap,
420 md::mdspan<const scalar_value_t<T>,
421 md::extents<std::size_t, md::dynamic_extent, 3>>
422 x,
423 FEkernel<T> auto kernel,
424 md::mdspan<const std::int32_t,
425 md::extents<std::size_t, md::dynamic_extent, 2, 2>>
426 facets,
427 std::tuple<mdspan2_t, int,
428 md::mdspan<const std::int32_t,
429 md::extents<std::size_t, md::dynamic_extent, 2, 2>>>
430 dofmap0,
431 fem::DofTransformKernel<T> auto P0,
432 std::tuple<mdspan2_t, int,
433 md::mdspan<const std::int32_t,
434 md::extents<std::size_t, md::dynamic_extent, 2, 2>>>
435 dofmap1,
436 fem::DofTransformKernel<T> auto P1T, std::span<const T> constants,
437 md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
438 md::dynamic_extent>>
439 coeffs,
440 std::span<const std::uint32_t> cell_info0,
441 std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
442 std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
443 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
444{
445 if (facets.empty())
446 return;
447
448 const auto [dmap0, bs0, facets0] = dofmap0;
449 const auto [dmap1, bs1, facets1] = dofmap1;
450
451 // Data structures used in assembly
452 using X = scalar_value_t<T>;
453 std::vector<X> cdofs(2 * x_dofmap.extent(1) * 3);
454 std::span<X> cdofs0(cdofs.data(), x_dofmap.extent(1) * 3);
455 std::span<X> cdofs1(cdofs.data() + x_dofmap.extent(1) * 3,
456 x_dofmap.extent(1) * 3);
457 std::vector<T> Ae, be;
458
459 // Temporaries for joint dofmaps
460 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
461
462 const int num_dofs0 = dmap0.extent(1);
463 const int num_dofs1 = dmap1.extent(1);
464 assert(facets0.size() == facets.size());
465 assert(facets1.size() == facets.size());
466 for (std::size_t f = 0; f < facets.extent(0); ++f)
467 {
468 // Cells in integration domain, test function domain and trial
469 // function domain meshes
470 std::array cells{facets(f, 0, 0), facets(f, 1, 0)};
471 std::array cells0{facets0(f, 0, 0), facets0(f, 1, 0)};
472 std::array cells1{facets1(f, 0, 0), facets1(f, 1, 0)};
473
474 // Local facet indices
475 std::array local_facet = {facets(f, 0, 1), facets(f, 1, 1)};
476
477 // Get cell geometry
478 auto x_dofs0 = md::submdspan(x_dofmap, cells[0], md::full_extent);
479 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
480 std::copy_n(&x(x_dofs0[i], 0), 3, std::next(cdofs0.begin(), 3 * i));
481 auto x_dofs1 = md::submdspan(x_dofmap, cells[1], md::full_extent);
482 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
483 std::copy_n(&x(x_dofs1[i], 0), 3, std::next(cdofs1.begin(), 3 * i));
484
485 // Get dof maps for cells and pack
486 std::span dmap0_cell0(dmap0.data_handle() + cells0[0] * num_dofs0,
487 num_dofs0);
488 std::span dmap0_cell1(dmap0.data_handle() + cells0[1] * num_dofs0,
489 num_dofs0);
490
491 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
492 std::ranges::copy(dmap0_cell0, dmapjoint0.begin());
493 std::ranges::copy(dmap0_cell1,
494 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
495
496 std::span dmap1_cell0(dmap1.data_handle() + cells1[0] * num_dofs1,
497 num_dofs1);
498 std::span dmap1_cell1(dmap1.data_handle() + cells1[1] * num_dofs1,
499 num_dofs1);
500
501 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
502 std::ranges::copy(dmap1_cell0, dmapjoint1.begin());
503 std::ranges::copy(dmap1_cell1,
504 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
505
506 // Check if bc is applied to cell0
507 bool has_bc = false;
508 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
509 {
510 for (int k = 0; k < bs1; ++k)
511 {
512 if (bc_markers1[bs1 * dmap1_cell0[j] + k])
513 {
514 has_bc = true;
515 break;
516 }
517 }
518 }
519
520 // Check if bc is applied to cell1
521 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
522 {
523 for (int k = 0; k < bs1; ++k)
524 {
525 if (bc_markers1[bs1 * dmap1_cell1[j] + k])
526 {
527 has_bc = true;
528 break;
529 }
530 }
531 }
532
533 if (!has_bc)
534 continue;
535
536 const int num_rows = bs0 * dmapjoint0.size();
537 const int num_cols = bs1 * dmapjoint1.size();
538
539 // Tabulate tensor
540 Ae.resize(num_rows * num_cols);
541 std::ranges::fill(Ae, 0);
542 std::array perm = perms.empty()
543 ? std::array<std::uint8_t, 2>{0, 0}
544 : std::array{perms(cells[0], local_facet[0]),
545 perms(cells[1], local_facet[1])};
546 kernel(Ae.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
547 local_facet.data(), perm.data(), nullptr);
548
549 std::span<T> _Ae(Ae);
550 std::span<T> sub_Ae0 = _Ae.subspan(bs0 * dmap0_cell0.size() * num_cols,
551 bs0 * dmap0_cell1.size() * num_cols);
552
553 P0(_Ae, cell_info0, cells0[0], num_cols);
554 P0(sub_Ae0, cell_info0, cells0[1], num_cols);
555 P1T(_Ae, cell_info1, cells1[0], num_rows);
556
557 for (int row = 0; row < num_rows; ++row)
558 {
559 // DOFs for dmap1 and cell1 are not stored contiguously in
560 // the block matrix, so each row needs a separate span access
561 std::span<T> sub_Ae1 = _Ae.subspan(
562 row * num_cols + bs1 * dmap1_cell0.size(), bs1 * dmap1_cell1.size());
563 P1T(sub_Ae1, cell_info1, cells1[1], 1);
564 }
565
566 be.resize(num_rows);
567 std::ranges::fill(be, 0);
568
569 // Compute b = b - A*b for cell0
570 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
571 {
572 for (int k = 0; k < bs1; ++k)
573 {
574 const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
575 if (bc_markers1[jj])
576 {
577 const T bc = bc_values1[jj];
578 const T _x0 = x0.empty() ? 0 : x0[jj];
579 // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
580 for (int m = 0; m < num_rows; ++m)
581 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
582 }
583 }
584 }
585
586 // Compute b = b - A*b for cell1
587 const int offset = bs1 * dmap1_cell0.size();
588 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
589 {
590 for (int k = 0; k < bs1; ++k)
591 {
592 const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
593 if (bc_markers1[jj])
594 {
595 const T bc = bc_values1[jj];
596 const T _x0 = x0.empty() ? 0 : x0[jj];
597 // be -= Ae.col(offset + bs1 * j + k) * alpha * (bc - x0[jj]);
598 for (int m = 0; m < num_rows; ++m)
599 {
600 be[m]
601 -= Ae[m * num_cols + offset + bs1 * j + k] * alpha * (bc - _x0);
602 }
603 }
604 }
605 }
606
607 for (std::size_t i = 0; i < dmap0_cell0.size(); ++i)
608 for (int k = 0; k < bs0; ++k)
609 b[bs0 * dmap0_cell0[i] + k] += be[bs0 * i + k];
610
611 const int offset_be = bs0 * dmap0_cell0.size();
612 for (std::size_t i = 0; i < dmap0_cell1.size(); ++i)
613 for (int k = 0; k < bs0; ++k)
614 b[bs0 * dmap0_cell1[i] + k] += be[offset_be + bs0 * i + k];
615 }
616}
617
641template <dolfinx::scalar T, int _bs = -1>
642void assemble_cells(
643 fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
644 md::mdspan<const scalar_value_t<T>,
645 md::extents<std::size_t, md::dynamic_extent, 3>>
646 x,
647 std::span<const std::int32_t> cells,
648 std::tuple<mdspan2_t, int, std::span<const std::int32_t>> dofmap,
649 FEkernel<T> auto kernel, std::span<const T> constants,
650 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
651 std::span<const std::uint32_t> cell_info0)
652{
653 if (cells.empty())
654 return;
655
656 const auto [dmap, bs, cells0] = dofmap;
657 assert(_bs < 0 or _bs == bs);
658
659 // Create data structures used in assembly
660 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
661 std::vector<T> be(bs * dmap.extent(1));
662 std::span<T> _be(be);
663
664 // Iterate over active cells
665 for (std::size_t index = 0; index < cells.size(); ++index)
666 {
667 // Integration domain celland test function cell
668 std::int32_t c = cells[index];
669 std::int32_t c0 = cells0[index];
670
671 // Get cell coordinates/geometry
672 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
673 for (std::size_t i = 0; i < x_dofs.size(); ++i)
674 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
675
676 // Tabulate vector for cell
677 std::ranges::fill(be, 0);
678 kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
679 nullptr, nullptr, nullptr);
680 P0(_be, cell_info0, c0, 1);
681
682 // Scatter cell vector to 'global' vector array
683 auto dofs = md::submdspan(dmap, c0, md::full_extent);
684 if constexpr (_bs > 0)
685 {
686 for (std::size_t i = 0; i < dofs.size(); ++i)
687 for (int k = 0; k < _bs; ++k)
688 b[_bs * dofs[i] + k] += be[_bs * i + k];
689 }
690 else
691 {
692 for (std::size_t i = 0; i < dofs.size(); ++i)
693 for (int k = 0; k < bs; ++k)
694 b[bs * dofs[i] + k] += be[bs * i + k];
695 }
696 }
697}
698
723template <dolfinx::scalar T, int _bs = -1>
724void assemble_exterior_facets(
725 fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
726 md::mdspan<const scalar_value_t<T>,
727 md::extents<std::size_t, md::dynamic_extent, 3>>
728 x,
729 md::mdspan<const std::int32_t,
730 std::extents<std::size_t, md::dynamic_extent, 2>>
731 facets,
732 std::tuple<mdspan2_t, int,
733 md::mdspan<const std::int32_t,
734 std::extents<std::size_t, md::dynamic_extent, 2>>>
735 dofmap,
736 FEkernel<T> auto fn, std::span<const T> constants,
737 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
738 std::span<const std::uint32_t> cell_info0,
739 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
740{
741 if (facets.empty())
742 return;
743
744 const auto [dmap, bs, facets0] = dofmap;
745 assert(_bs < 0 or _bs == bs);
746
747 // Create data structures used in assembly
748 const int num_dofs = dmap.extent(1);
749 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
750 std::vector<T> be(bs * num_dofs);
751 std::span<T> _be(be);
752 assert(facets0.size() == facets.size());
753 for (std::size_t f = 0; f < facets.extent(0); ++f)
754 {
755 // Cell in the integration domain, local facet index relative to the
756 // integration domain cell, and cell in the test function mesh
757 std::int32_t cell = facets(f, 0);
758 std::int32_t local_facet = facets(f, 1);
759 std::int32_t cell0 = facets0(f, 0);
760
761 // Get cell coordinates/geometry
762 auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
763 for (std::size_t i = 0; i < x_dofs.size(); ++i)
764 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
765
766 // Permutations
767 std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
768
769 // Tabulate element vector
770 std::ranges::fill(be, 0);
771 fn(be.data(), &coeffs(f, 0), constants.data(), cdofs.data(), &local_facet,
772 &perm, nullptr);
773
774 P0(_be, cell_info0, cell0, 1);
775
776 // Add element vector to global vector
777 auto dofs = md::submdspan(dmap, cell0, md::full_extent);
778 if constexpr (_bs > 0)
779 {
780 for (std::size_t i = 0; i < dofs.size(); ++i)
781 for (int k = 0; k < _bs; ++k)
782 b[_bs * dofs[i] + k] += be[_bs * i + k];
783 }
784 else
785 {
786 for (std::size_t i = 0; i < dofs.size(); ++i)
787 for (int k = 0; k < bs; ++k)
788 b[bs * dofs[i] + k] += be[bs * i + k];
789 }
790 }
791}
792
816template <dolfinx::scalar T, int _bs = -1>
817void assemble_interior_facets(
818 fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
819 md::mdspan<const scalar_value_t<T>,
820 md::extents<std::size_t, md::dynamic_extent, 3>>
821 x,
822 md::mdspan<const std::int32_t,
823 std::extents<std::size_t, md::dynamic_extent, 2, 2>>
824 facets,
825 std::tuple<const DofMap&, int,
826 md::mdspan<const std::int32_t,
827 std::extents<std::size_t, md::dynamic_extent, 2, 2>>>
828 dofmap,
829 FEkernel<T> auto fn, std::span<const T> constants,
830 md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
831 md::dynamic_extent>>
832 coeffs,
833 std::span<const std::uint32_t> cell_info0,
834 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
835{
836 if (facets.empty())
837 return;
838
839 const auto [dmap, bs, facets0] = dofmap;
840 assert(_bs < 0 or _bs == bs);
841
842 // Create data structures used in assembly
843 using X = scalar_value_t<T>;
844 std::vector<X> cdofs(2 * x_dofmap.extent(1) * 3);
845 std::span<X> cdofs0(cdofs.data(), x_dofmap.extent(1) * 3);
846 std::span<X> cdofs1(cdofs.data() + x_dofmap.extent(1) * 3,
847 x_dofmap.extent(1) * 3);
848 std::vector<T> be;
849
850 assert(facets0.size() == facets.size());
851 for (std::size_t f = 0; f < facets.extent(0); ++f)
852 {
853 // Cells in integration domain and test function domain meshes
854 std::array<std::int32_t, 2> cells{facets(f, 0, 0), facets(f, 1, 0)};
855 std::array<std::int32_t, 2> cells0{facets0(f, 0, 0), facets0(f, 1, 0)};
856
857 // Local facet indices
858 std::array<std::int32_t, 2> local_facet{facets(f, 0, 1), facets(f, 1, 1)};
859
860 // Get cell geometry
861 auto x_dofs0 = md::submdspan(x_dofmap, cells[0], md::full_extent);
862 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
863 std::copy_n(&x(x_dofs0[i], 0), 3, std::next(cdofs0.begin(), 3 * i));
864 auto x_dofs1 = md::submdspan(x_dofmap, cells[1], md::full_extent);
865 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
866 std::copy_n(&x(x_dofs1[i], 0), 3, std::next(cdofs1.begin(), 3 * i));
867
868 // Get dofmaps for cells
869 std::span dmap0 = dmap.cell_dofs(cells0[0]);
870 std::span 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 = perms.empty()
876 ? std::array<std::uint8_t, 2>{0, 0}
877 : std::array{perms(cells[0], local_facet[0]),
878 perms(cells[1], local_facet[1])};
879 fn(be.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
880 local_facet.data(), perm.data(), nullptr);
881
882 std::span<T> _be(be);
883 std::span<T> sub_be = _be.subspan(bs * dmap0.size(), bs * dmap1.size());
884
885 P0(be, cell_info0, cells0[0], 1);
886 P0(sub_be, cell_info0, cells0[1], 1);
887
888 // Add element vector to global vector
889 if constexpr (_bs > 0)
890 {
891 for (std::size_t i = 0; i < dmap0.size(); ++i)
892 for (int k = 0; k < _bs; ++k)
893 b[_bs * dmap0[i] + k] += be[_bs * i + k];
894 for (std::size_t i = 0; i < dmap1.size(); ++i)
895 for (int k = 0; k < _bs; ++k)
896 b[_bs * dmap1[i] + k] += be[_bs * (i + dmap0.size()) + k];
897 }
898 else
899 {
900 for (std::size_t i = 0; i < dmap0.size(); ++i)
901 for (int k = 0; k < bs; ++k)
902 b[bs * dmap0[i] + k] += be[bs * i + k];
903 for (std::size_t i = 0; i < dmap1.size(); ++i)
904 for (int k = 0; k < bs; ++k)
905 b[bs * dmap1[i] + k] += be[bs * (i + dmap0.size()) + k];
906 }
907 }
908}
909
926template <dolfinx::scalar T, std::floating_point U>
927void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
928 md::mdspan<const scalar_value_t<T>,
929 md::extents<std::size_t, md::dynamic_extent, 3>>
930 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
942 // Test function mesh
943 auto mesh0 = a.function_spaces().at(0)->mesh();
944 assert(mesh0);
945
946 // Trial function mesh
947 auto mesh1 = a.function_spaces().at(1)->mesh();
948 assert(mesh1);
949
950 // Get dofmap for columns and rows of a
951 assert(a.function_spaces().at(0));
952 assert(a.function_spaces().at(1));
953 auto dofmap0 = a.function_spaces()[0]->dofmap()->map();
954 const int bs0 = a.function_spaces()[0]->dofmap()->bs();
955 auto element0 = a.function_spaces()[0]->element();
956 auto dofmap1 = a.function_spaces()[1]->dofmap()->map();
957 const int bs1 = a.function_spaces()[1]->dofmap()->bs();
958 auto element1 = a.function_spaces()[1]->element();
959 assert(element0);
960
961 std::span<const std::uint32_t> cell_info0;
962 std::span<const std::uint32_t> cell_info1;
963 // TODO: Check for each element instead
964 if (element0->needs_dof_transformations()
965 or element1->needs_dof_transformations() or a.needs_facet_permutations())
966 {
967 mesh0->topology_mutable()->create_entity_permutations();
968 mesh1->topology_mutable()->create_entity_permutations();
969 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
970 cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
971 }
972
973 fem::DofTransformKernel<T> auto P0
974 = element0->template dof_transformation_fn<T>(doftransform::standard);
975 fem::DofTransformKernel<T> auto P1T
976 = element1->template dof_transformation_right_fn<T>(
978
979 for (int i : a.integral_ids(IntegralType::cell))
980 {
981 auto kernel = a.kernel(IntegralType::cell, i, 0);
982 assert(kernel);
983 auto& [_coeffs, cstride] = coefficients.at({IntegralType::cell, i});
984 std::span cells = a.domain(IntegralType::cell, i, 0);
985 std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, 0);
986 std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, 0);
987 assert(_coeffs.size() == cells.size() * cstride);
988 auto coeffs = md::mdspan(_coeffs.data(), cells.size(), cstride);
989 if (bs0 == 1 and bs1 == 1)
990 {
991 _lift_bc_cells<T, 1, 1>(
992 b, x_dofmap, x, kernel, cells, {dofmap0, bs0, cells0}, P0,
993 {dofmap1, bs1, cells1}, P1T, constants, coeffs, cell_info0,
994 cell_info1, bc_values1, bc_markers1, x0, alpha);
995 }
996 else if (bs0 == 3 and bs1 == 3)
997 {
998 _lift_bc_cells<T, 3, 3>(
999 b, x_dofmap, x, kernel, cells, {dofmap0, bs0, cells0}, P0,
1000 {dofmap1, bs1, cells1}, P1T, constants, coeffs, cell_info0,
1001 cell_info1, bc_values1, bc_markers1, x0, alpha);
1002 }
1003 else
1004 {
1005 _lift_bc_cells(b, x_dofmap, x, kernel, cells, {dofmap0, bs0, cells0}, P0,
1006 {dofmap1, bs1, cells1}, P1T, constants, coeffs, cell_info0,
1007 cell_info1, bc_values1, bc_markers1, x0, alpha);
1008 }
1009 }
1010
1011 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
1012 if (a.needs_facet_permutations())
1013 {
1014 mesh::CellType cell_type = mesh->topology()->cell_type();
1015 int num_facets_per_cell
1016 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
1017 mesh->topology_mutable()->create_entity_permutations();
1018 const std::vector<std::uint8_t>& p
1019 = mesh->topology()->get_facet_permutations();
1020 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
1021 num_facets_per_cell);
1022 }
1023
1024 for (int i : a.integral_ids(IntegralType::exterior_facet))
1025 {
1026 auto kernel = a.kernel(IntegralType::exterior_facet, i, 0);
1027 assert(kernel);
1028 auto& [coeffs, cstride]
1029 = coefficients.at({IntegralType::exterior_facet, i});
1030
1031 using mdspanx2_t
1032 = md::mdspan<const std::int32_t,
1033 md::extents<std::size_t, md::dynamic_extent, 2>>;
1034 std::span f = a.domain(IntegralType::exterior_facet, i, 0);
1035 mdspanx2_t facets(f.data(), f.size() / 2, 2);
1036 std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0);
1037 mdspanx2_t facets0(f0.data(), f0.size() / 2, 2);
1038 std::span f1 = a.domain_arg(IntegralType::exterior_facet, 1, i, 0);
1039 mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
1040 assert(coeffs.size() == facets.extent(0) * cstride);
1041 _lift_bc_exterior_facets(
1042 b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
1043 {dofmap1, bs1, facets1}, P1T, constants,
1044 md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
1045 cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
1046 }
1047
1048 for (int i : a.integral_ids(IntegralType::interior_facet))
1049 {
1050 auto kernel = a.kernel(IntegralType::interior_facet, i, 0);
1051 assert(kernel);
1052 auto& [coeffs, cstride]
1053 = coefficients.at({IntegralType::interior_facet, i});
1054
1055 using mdspanx22_t
1056 = md::mdspan<const std::int32_t,
1057 md::extents<std::size_t, md::dynamic_extent, 2, 2>>;
1058 using mdspanx2x_t
1059 = md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
1060 md::dynamic_extent>>;
1061 std::span f = a.domain(IntegralType::interior_facet, i, 0);
1062 mdspanx22_t facets(f.data(), f.size() / 4, 2, 2);
1063 std::span f0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0);
1064 mdspanx22_t facets0(f0.data(), f0.size() / 4, 2, 2);
1065 std::span f1 = a.domain_arg(IntegralType::interior_facet, 1, i, 0);
1066 mdspanx22_t facets1(f1.data(), f1.size() / 4, 2, 2);
1067 _lift_bc_interior_facets(
1068 b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
1069 {dofmap1, bs1, facets1}, P1T, constants,
1070 mdspanx2x_t(coeffs.data(), facets.extent(0), 2, cstride), cell_info0,
1071 cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
1072 }
1073}
1074
1095template <dolfinx::scalar T, std::floating_point U>
1096void apply_lifting(
1097 std::span<T> b,
1098 std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> a,
1099 const std::vector<std::span<const T>>& constants,
1100 const std::vector<std::map<std::pair<IntegralType, int>,
1101 std::pair<std::span<const T>, int>>>& coeffs,
1102 const std::vector<
1103 std::vector<std::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
1104 const std::vector<std::span<const T>>& x0, T alpha)
1105{
1106 if (!x0.empty() and x0.size() != a.size())
1107 {
1108 throw std::runtime_error(
1109 "Mismatch in size between x0 and bilinear form in assembler.");
1110 }
1111
1112 if (a.size() != bcs1.size())
1113 {
1114 throw std::runtime_error(
1115 "Mismatch in size between a and bcs in assembler.");
1116 }
1117
1118 for (std::size_t j = 0; j < a.size(); ++j)
1119 {
1120 std::vector<std::int8_t> bc_markers1;
1121 std::vector<T> bc_values1;
1122 if (a[j] and !bcs1[j].empty())
1123 {
1124 // Extract data from mesh
1125 std::shared_ptr<const mesh::Mesh<U>> mesh = a[j]->get().mesh();
1126 if (!mesh)
1127 throw std::runtime_error("Unable to extract a mesh.");
1128 mdspan2_t x_dofmap = mesh->geometry().dofmap();
1129 std::span _x = mesh->geometry().x();
1130 md::mdspan<const scalar_value_t<T>,
1131 md::extents<std::size_t, md::dynamic_extent, 3>>
1132 x(_x.data(), _x.size() / 3, 3);
1133
1134 assert(a[j]->get().function_spaces().at(0));
1135 auto V1 = a[j]->get().function_spaces()[1];
1136 assert(V1);
1137 auto map1 = V1->dofmap()->index_map;
1138 const int bs1 = V1->dofmap()->index_map_bs();
1139 assert(map1);
1140 const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
1141 bc_markers1.assign(crange, false);
1142 bc_values1.assign(crange, 0);
1143 for (auto& bc : bcs1[j])
1144 {
1145 bc.get().mark_dofs(bc_markers1);
1146 bc.get().set(bc_values1, std::nullopt, 1);
1147 }
1148
1149 if (!x0.empty())
1150 {
1151 lift_bc<T>(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j],
1152 bc_values1, bc_markers1, x0[j], alpha);
1153 }
1154 else
1155 {
1156 lift_bc<T>(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j],
1157 bc_values1, bc_markers1, std::span<const T>(), alpha);
1158 }
1159 }
1160 }
1161}
1162
1170template <dolfinx::scalar T, std::floating_point U>
1171void assemble_vector(
1172 std::span<T> b, const Form<T, U>& L,
1173 md::mdspan<const scalar_value_t<T>,
1174 md::extents<std::size_t, md::dynamic_extent, 3>>
1175 x,
1176 std::span<const T> constants,
1177 const std::map<std::pair<IntegralType, int>,
1178 std::pair<std::span<const T>, int>>& coefficients)
1179{
1180 // Integration domain mesh
1181 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1182 assert(mesh);
1183
1184 // Test function mesh
1185 auto mesh0 = L.function_spaces().at(0)->mesh();
1186 assert(mesh0);
1187
1188 const int num_cell_types = mesh->topology()->cell_types().size();
1189 for (int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
1190 {
1191 // Geometry dofmap and data
1192 mdspan2_t x_dofmap = mesh->geometry().dofmap(cell_type_idx);
1193
1194 // Get dofmap data
1195 assert(L.function_spaces().at(0));
1196 auto element = L.function_spaces().at(0)->elements(cell_type_idx);
1197 assert(element);
1198 std::shared_ptr<const fem::DofMap> dofmap
1199 = L.function_spaces().at(0)->dofmaps(cell_type_idx);
1200 assert(dofmap);
1201 auto dofs = dofmap->map();
1202 const int bs = dofmap->bs();
1203
1204 fem::DofTransformKernel<T> auto P0
1205 = element->template dof_transformation_fn<T>(doftransform::standard);
1206
1207 std::span<const std::uint32_t> cell_info0;
1208 if (element->needs_dof_transformations() or L.needs_facet_permutations())
1209 {
1210 mesh0->topology_mutable()->create_entity_permutations();
1211 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
1212 }
1213
1214 for (int i : L.integral_ids(IntegralType::cell))
1215 {
1216 auto fn = L.kernel(IntegralType::cell, i, cell_type_idx);
1217 assert(fn);
1218 std::span cells = L.domain(IntegralType::cell, i, cell_type_idx);
1219 std::span cells0 = L.domain_arg(IntegralType::cell, 0, i, cell_type_idx);
1220 auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
1221 assert(cells.size() * cstride == coeffs.size());
1222 if (bs == 1)
1223 {
1224 impl::assemble_cells<T, 1>(
1225 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1226 md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
1227 }
1228 else if (bs == 3)
1229 {
1230 impl::assemble_cells<T, 3>(
1231 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1232 md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
1233 }
1234 else
1235 {
1236 impl::assemble_cells(
1237 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1238 md::mdspan(coeffs.data(), cells.size(), cstride), cell_info0);
1239 }
1240 }
1241
1242 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
1243 if (L.needs_facet_permutations())
1244 {
1245 mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx];
1246 int num_facets_per_cell
1247 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
1248 mesh->topology_mutable()->create_entity_permutations();
1249 const std::vector<std::uint8_t>& p
1250 = mesh->topology()->get_facet_permutations();
1251 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
1252 num_facets_per_cell);
1253 }
1254
1255 using mdspanx2_t
1256 = md::mdspan<const std::int32_t,
1257 md::extents<std::size_t, md::dynamic_extent, 2>>;
1258
1259 for (int i : L.integral_ids(IntegralType::exterior_facet))
1260 {
1261 auto fn = L.kernel(IntegralType::exterior_facet, i, 0);
1262 assert(fn);
1263 auto& [coeffs, cstride]
1264 = coefficients.at({IntegralType::exterior_facet, i});
1265 std::span f = L.domain(IntegralType::exterior_facet, i, 0);
1266 mdspanx2_t facets(f.data(), f.size() / 2, 2);
1267 std::span f1 = L.domain_arg(IntegralType::exterior_facet, 0, i, 0);
1268 mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
1269 assert((facets.size() / 2) * cstride == coeffs.size());
1270 if (bs == 1)
1271 {
1272 impl::assemble_exterior_facets<T, 1>(
1273 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1274 md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
1275 perms);
1276 }
1277 else if (bs == 3)
1278 {
1279 impl::assemble_exterior_facets<T, 3>(
1280 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1281 md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
1282 perms);
1283 }
1284 else
1285 {
1286 impl::assemble_exterior_facets(
1287 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1288 md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
1289 perms);
1290 }
1291 }
1292
1293 for (int i : L.integral_ids(IntegralType::interior_facet))
1294 {
1295 using mdspanx22_t
1296 = md::mdspan<const std::int32_t,
1297 md::extents<std::size_t, md::dynamic_extent, 2, 2>>;
1298 using mdspanx2x_t
1299 = md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
1300 md::dynamic_extent>>;
1301
1302 auto fn = L.kernel(IntegralType::interior_facet, i, 0);
1303 assert(fn);
1304 auto& [coeffs, cstride]
1305 = coefficients.at({IntegralType::interior_facet, i});
1306 std::span facets = L.domain(IntegralType::interior_facet, i, 0);
1307 std::span facets1 = L.domain_arg(IntegralType::interior_facet, 0, i, 0);
1308 assert((facets.size() / 4) * 2 * cstride == coeffs.size());
1309 if (bs == 1)
1310 {
1311 impl::assemble_interior_facets<T, 1>(
1312 P0, b, x_dofmap, x,
1313 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1314 {*dofmap, bs,
1315 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1316 fn, constants,
1317 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1318 cell_info0, perms);
1319 }
1320 else if (bs == 3)
1321 {
1322 impl::assemble_interior_facets<T, 3>(
1323 P0, b, x_dofmap, x,
1324 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1325 {*dofmap, bs,
1326 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1327 fn, constants,
1328 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1329 cell_info0, perms);
1330 }
1331 else
1332 {
1333 impl::assemble_interior_facets(
1334 P0, b, x_dofmap, x,
1335 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1336 {*dofmap, bs,
1337 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1338 fn, constants,
1339 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1340 cell_info0, perms);
1341 }
1342 }
1343 }
1344}
1345
1352template <dolfinx::scalar T, std::floating_point U>
1353void assemble_vector(
1354 std::span<T> b, const Form<T, U>& L, std::span<const T> constants,
1355 const std::map<std::pair<IntegralType, int>,
1356 std::pair<std::span<const T>, int>>& coefficients)
1357{
1358 using mdspanx3_t
1359 = md::mdspan<const scalar_value_t<T>,
1360 md::extents<std::size_t, md::dynamic_extent, 3>>;
1361
1362 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1363 assert(mesh);
1364 auto x = mesh->geometry().x();
1365 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
1366 {
1367 assemble_vector(b, L, mdspanx3_t(x.data(), x.size() / 3, 3), constants,
1368 coefficients);
1369 }
1370 else
1371 {
1372 std::vector<scalar_value_t<T>> _x(x.begin(), x.end());
1373 assemble_vector(b, L, mdspanx3_t(_x.data(), _x.size() / 3, 3), constants,
1374 coefficients);
1375 }
1376}
1377} // 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
@ interior_facet
Interior facet.
Definition Form.h:39
@ cell
Cell.
Definition Form.h:37
@ exterior_facet
Exterior facet.
Definition Form.h:38
int cell_num_entities(CellType type, int dim)
Definition cell_types.cpp:139
CellType
Cell type identifier.
Definition cell_types.h:22