10#include "DirichletBC.h"
16#include <basix/mdspan.hpp>
18#include <dolfinx/common/IndexMap.h>
19#include <dolfinx/mesh/Geometry.h>
20#include <dolfinx/mesh/Mesh.h>
21#include <dolfinx/mesh/Topology.h>
30template <dolfinx::scalar T, std::
floating_po
int U>
34namespace dolfinx::fem::impl
37using mdspan2_t = md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>;
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>
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>>
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)
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);
105 const int num_rows = bs0 * dmap0.extent(1);
106 const int num_cols = bs1 * dmap1.extent(1);
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)
117 std::int32_t c =
cells[index];
118 std::int32_t c0 = cells0[index];
119 std::int32_t c1 = cells1[index];
122 auto dofs1 = md::submdspan(dmap1, c1, md::full_extent);
126 for (std::size_t j = 0; j < dofs1.size(); ++j)
128 if constexpr (_bs1 > 0)
130 for (
int k = 0; k < _bs1; ++k)
132 assert(_bs1 * dofs1[j] + k < (
int)bc_markers1.size());
133 if (bc_markers1[_bs1 * dofs1[j] + k])
142 for (
int k = 0; k < bs1; ++k)
144 assert(bs1 * dofs1[j] + k < (
int)bc_markers1.size());
145 if (bc_markers1[bs1 * dofs1[j] + k])
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));
163 auto dofs0 = md::submdspan(dmap0, c0, md::full_extent);
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);
172 std::ranges::fill(be, 0);
173 for (std::size_t j = 0; j < dofs1.size(); ++j)
175 if constexpr (_bs1 > 0)
177 for (
int k = 0; k < _bs1; ++k)
179 const std::int32_t jj = _bs1 * dofs1[j] + k;
180 assert(jj < (
int)bc_markers1.size());
183 const T bc = bc_values1[jj];
184 const T _x0 = x0.empty() ? 0 : x0[jj];
187 for (
int m = 0; m < num_rows; ++m)
188 be[m] -= Ae[m * num_cols + _bs1 * j + k] * alpha * (bc - _x0);
194 for (
int k = 0; k < bs1; ++k)
196 const std::int32_t jj = bs1 * dofs1[j] + k;
197 assert(jj < (
int)bc_markers1.size());
200 const T bc = bc_values1[jj];
201 const T _x0 = x0.empty() ? 0 : x0[jj];
203 for (
int m = 0; m < num_rows; ++m)
204 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
210 for (std::size_t i = 0; i < dofs0.size(); ++i)
212 if constexpr (_bs0 > 0)
214 for (
int k = 0; k < _bs0; ++k)
215 b[_bs0 * dofs0[i] + k] += be[_bs0 * i + k];
219 for (
int k = 0; k < bs0; ++k)
220 b[bs0 * dofs0[i] + k] += be[bs0 * i + k];
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>>
271 FEkernel<T>
auto kernel,
272 md::mdspan<
const std::int32_t,
273 md::extents<std::size_t, md::dynamic_extent, 2>>
275 std::tuple<mdspan2_t,
int,
276 md::mdspan<
const std::int32_t,
277 md::extents<std::size_t, md::dynamic_extent, 2>>>
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>>>
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)
294 const auto [dmap0, bs0, facets0] = dofmap0;
295 const auto [dmap1, bs1, facets1] = dofmap1;
297 const int num_rows = bs0 * dmap0.extent(1);
298 const int num_cols = bs1 * dmap1.extent(1);
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)
309 std::int32_t
cell = facets(index, 0);
310 std::int32_t cell0 = facets0(index, 0);
311 std::int32_t cell1 = facets1(index, 0);
314 std::int32_t local_facet = facets(index, 1);
317 auto dofs1 = md::submdspan(dmap1, cell1, md::full_extent);
321 for (std::size_t j = 0; j < dofs1.size(); ++j)
323 for (
int k = 0; k < bs1; ++k)
325 if (bc_markers1[bs1 * dofs1[j] + k])
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));
342 auto dofs0 = md::submdspan(dmap0, cell0, md::full_extent);
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);
353 std::ranges::fill(be, 0);
354 for (std::size_t j = 0; j < dofs1.size(); ++j)
356 for (
int k = 0; k < bs1; ++k)
358 const std::int32_t jj = bs1 * dofs1[j] + k;
361 const T bc = bc_values1[jj];
362 const T _x0 = x0.empty() ? 0 : x0[jj];
364 for (
int m = 0; m < num_rows; ++m)
365 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
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];
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>>
426 FEkernel<T>
auto kernel,
427 md::mdspan<
const std::int32_t,
428 md::extents<std::size_t, md::dynamic_extent, 2, 2>>
430 std::tuple<mdspan2_t,
int,
431 md::mdspan<
const std::int32_t,
432 md::extents<std::size_t, md::dynamic_extent, 2, 2>>>
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>>>
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,
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)
451 const auto [dmap0, bs0, facets0] = dofmap0;
452 const auto [dmap1, bs1, facets1] = dofmap1;
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;
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);
468 std::vector<std::int32_t> dmapjoint0(2 * num_dofs0);
469 std::vector<std::int32_t> dmapjoint1(2 * num_dofs1);
471 assert(facets0.size() == facets.size());
472 assert(facets1.size() == facets.size());
473 for (std::size_t f = 0; f < facets.extent(0); ++f)
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)};
482 std::array local_facet = {facets(f, 0, 1), facets(f, 1, 1)};
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));
496 std::span dmap0_cell0
498 ? std::span(dmap0.data_handle() + cells0[0] * num_dofs0,
500 : std::span<const std::int32_t>();
501 std::span dmap0_cell1
503 ? std::span(dmap0.data_handle() + cells0[1] * num_dofs0,
505 : std::span<const std::int32_t>();
507 std::ranges::copy(dmap0_cell0, dmapjoint0.begin());
508 std::ranges::copy(dmap0_cell1, std::next(dmapjoint0.begin(), num_dofs0));
511 std::span<const std::int32_t> dmap1_cell0
513 ? std::span(dmap1.data_handle() + cells1[0] * num_dofs1,
515 : std::span<const std::int32_t>();
516 std::span<const std::int32_t> dmap1_cell1
518 ? std::span(dmap1.data_handle() + cells1[1] * num_dofs1,
520 : std::span<const std::int32_t>();
522 std::ranges::copy(dmap1_cell0, dmapjoint1.begin());
523 std::ranges::copy(dmap1_cell1, std::next(dmapjoint1.begin(), num_dofs1));
527 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
529 for (
int k = 0; k < bs1; ++k)
531 if (bc_markers1[bs1 * dmap1_cell0[j] + k])
540 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
542 for (
int k = 0; k < bs1; ++k)
544 if (bc_markers1[bs1 * dmap1_cell1[j] + k])
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);
565 P0(Ae, cell_info0, cells0[0], num_cols);
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);
573 P1T(Ae, cell_info1, cells1[0], num_rows);
577 for (
int row = 0; row < num_rows; ++row)
581 std::span sub_Ae1(Ae.data() + row * num_cols + bs1 * num_dofs1,
583 P1T(sub_Ae1, cell_info1, cells1[1], 1);
587 std::ranges::fill(be, 0);
590 for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
592 for (
int k = 0; k < bs1; ++k)
594 const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
597 const T bc = bc_values1[jj];
598 const T _x0 = x0.empty() ? 0 : x0[jj];
600 for (
int m = 0; m < num_rows; ++m)
601 be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
607 const int offset = bs1 * num_dofs1;
608 for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
610 for (
int k = 0; k < bs1; ++k)
612 const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
615 const T bc = bc_values1[jj];
616 const T _x0 = x0.empty() ? 0 : x0[jj];
618 for (
int m = 0; m < num_rows; ++m)
621 -= Ae[m * num_cols + offset + bs1 * j + k] * alpha * (bc - _x0);
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];
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];
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>
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>>
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)
679 const auto [dmap, bs, cells0] = dofmap;
680 assert(_bs < 0 or _bs == bs);
683 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
684 std::vector<T> be(bs * dmap.extent(1));
687 for (std::size_t index = 0; index <
cells.size(); ++index)
690 std::int32_t c =
cells[index];
691 std::int32_t c0 = cells0[index];
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));
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);
705 auto dofs = md::submdspan(dmap, c0, md::full_extent);
706 if constexpr (_bs > 0)
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];
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];
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>>
753 md::mdspan<
const std::int32_t,
754 std::extents<std::size_t, md::dynamic_extent, 2>>
756 std::tuple<mdspan2_t,
int,
757 md::mdspan<
const std::int32_t,
758 std::extents<std::size_t, md::dynamic_extent, 2>>>
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)
768 const auto [dmap, bs, facets0] = dofmap;
769 assert(_bs < 0 or _bs == bs);
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)
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);
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));
790 std::uint8_t perm = perms.empty() ? 0 : perms(
cell, local_facet);
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);
799 auto dofs = md::submdspan(dmap, cell0, md::full_extent);
800 if constexpr (_bs > 0)
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];
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];
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>>
848 md::mdspan<
const std::int32_t,
849 std::extents<std::size_t, md::dynamic_extent, 2, 2>>
851 std::tuple<
const DofMap&,
int,
852 md::mdspan<
const std::int32_t,
853 std::extents<std::size_t, md::dynamic_extent, 2, 2>>>
855 FEkernel<T>
auto kernel, std::span<const T> constants,
856 md::mdspan<
const T, md::extents<std::size_t, md::dynamic_extent, 2,
859 std::span<const std::uint32_t> cell_info0,
860 md::mdspan<
const std::uint8_t, md::dextents<std::size_t, 2>> perms)
862 using X = scalar_value_t<T>;
867 const auto [dmap, bs, facets0] = dofmap;
868 assert(_bs < 0 or _bs == bs);
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);
876 const std::size_t dmap_size = dmap.map().extent(1);
877 std::vector<T> be(bs * 2 * dmap_size);
879 assert(facets0.size() == facets.size());
880 for (std::size_t f = 0; f < facets.extent(0); ++f)
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)};
887 std::array<std::int32_t, 2> local_facet{facets(f, 0, 1), facets(f, 1, 1)};
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));
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>();
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);
915 P0(be, cell_info0, cells0[0], 1);
918 std::span sub_be(be.data() + bs * dmap_size, bs * dmap_size);
919 P0(sub_be, cell_info0, cells0[1], 1);
923 if constexpr (_bs > 0)
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];
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];
960template <
typename V, std::floating_point U,
961 dolfinx::scalar T =
typename std::remove_cvref_t<V>::value_type>
962 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
963void lift_bc(V&& b,
const Form<T, U>& a, mdspan2_t x_dofmap,
964 md::mdspan<
const scalar_value_t<T>,
965 md::extents<std::size_t, md::dynamic_extent, 3>>
967 std::span<const T> constants,
968 const std::map<std::pair<IntegralType, int>,
969 std::pair<std::span<const T>,
int>>& coefficients,
970 std::span<const T> bc_values1,
971 std::span<const std::int8_t> bc_markers1, std::span<const T> x0,
975 std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
979 auto mesh0 = a.function_spaces().at(0)->mesh();
983 auto mesh1 = a.function_spaces().at(1)->mesh();
987 assert(a.function_spaces().at(0));
988 assert(a.function_spaces().at(1));
989 auto dofmap0 = a.function_spaces()[0]->dofmap()->map();
990 const int bs0 = a.function_spaces()[0]->dofmap()->bs();
991 auto element0 = a.function_spaces()[0]->element();
992 auto dofmap1 = a.function_spaces()[1]->dofmap()->map();
993 const int bs1 = a.function_spaces()[1]->dofmap()->bs();
994 auto element1 = a.function_spaces()[1]->element();
997 std::span<const std::uint32_t> cell_info0;
998 std::span<const std::uint32_t> cell_info1;
1000 if (element0->needs_dof_transformations()
1001 or element1->needs_dof_transformations() or a.needs_facet_permutations())
1003 mesh0->topology_mutable()->create_entity_permutations();
1004 mesh1->topology_mutable()->create_entity_permutations();
1005 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
1006 cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
1009 fem::DofTransformKernel<T>
auto P0
1011 fem::DofTransformKernel<T>
auto P1T
1012 = element1->template dof_transformation_right_fn<T>(
1023 assert(_coeffs.size() ==
cells.size() * cstride);
1024 auto coeffs = md::mdspan(_coeffs.data(),
cells.size(), cstride);
1025 if (bs0 == 1 and bs1 == 1)
1027 _lift_bc_cells<1, 1>(b, x_dofmap, x, kernel, cells,
1028 {dofmap0, bs0, cells0}, P0, {dofmap1, bs1, cells1},
1029 P1T, constants, coeffs, cell_info0, cell_info1,
1030 bc_values1, bc_markers1, x0, alpha);
1032 else if (bs0 == 3 and bs1 == 3)
1034 _lift_bc_cells<3, 3>(b, x_dofmap, x, kernel, cells,
1035 {dofmap0, bs0, cells0}, P0, {dofmap1, bs1, cells1},
1036 P1T, constants, coeffs, cell_info0, cell_info1,
1037 bc_values1, bc_markers1, x0, alpha);
1041 _lift_bc_cells(b, x_dofmap, x, kernel, cells, {dofmap0, bs0, cells0}, P0,
1042 {dofmap1, bs1, cells1}, P1T, constants, coeffs, cell_info0,
1043 cell_info1, bc_values1, bc_markers1, x0, alpha);
1047 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
1048 if (a.needs_facet_permutations())
1051 int num_facets_per_cell
1053 mesh->topology_mutable()->create_entity_permutations();
1054 const std::vector<std::uint8_t>& p
1055 = mesh->topology()->get_facet_permutations();
1056 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
1057 num_facets_per_cell);
1064 auto& [coeffs, cstride]
1068 = md::mdspan<
const std::int32_t,
1069 md::extents<std::size_t, md::dynamic_extent, 2>>;
1071 mdspanx2_t facets(f.data(), f.size() / 2, 2);
1073 mdspanx2_t facets0(f0.data(), f0.size() / 2, 2);
1075 mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
1076 assert(coeffs.size() == facets.extent(0) * cstride);
1077 _lift_bc_exterior_facets(
1078 b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
1079 {dofmap1, bs1, facets1}, P1T, constants,
1080 md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
1081 cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
1088 auto& [coeffs, cstride]
1092 = md::mdspan<
const std::int32_t,
1093 md::extents<std::size_t, md::dynamic_extent, 2, 2>>;
1095 = md::mdspan<
const T, md::extents<std::size_t, md::dynamic_extent, 2,
1096 md::dynamic_extent>>;
1098 mdspanx22_t facets(f.data(), f.size() / 4, 2, 2);
1100 mdspanx22_t facets0(f0.data(), f0.size() / 4, 2, 2);
1102 mdspanx22_t facets1(f1.data(), f1.size() / 4, 2, 2);
1103 _lift_bc_interior_facets(
1104 b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0,
1105 {dofmap1, bs1, facets1}, P1T, constants,
1106 mdspanx2x_t(coeffs.data(), facets.extent(0), 2, cstride), cell_info0,
1107 cell_info1, bc_values1, bc_markers1, x0, alpha, perms);
1131template <
typename V, std::floating_point U,
1132 dolfinx::scalar T =
typename std::remove_cvref_t<V>::value_type>
1133 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
1136 std::vector<std::optional<std::reference_wrapper<
const Form<T, U>>>> a,
1137 const std::vector<std::span<const T>>& constants,
1138 const std::vector<std::map<std::pair<IntegralType, int>,
1139 std::pair<std::span<const T>,
int>>>& coeffs,
1141 std::vector<std::reference_wrapper<
const DirichletBC<T, U>>>>& bcs1,
1142 const std::vector<std::span<const T>>& x0, T alpha)
1144 if (!x0.empty() and x0.size() != a.size())
1146 throw std::runtime_error(
1147 "Mismatch in size between x0 and bilinear form in assembler.");
1150 if (a.size() != bcs1.size())
1152 throw std::runtime_error(
1153 "Mismatch in size between a and bcs in assembler.");
1156 for (std::size_t j = 0; j < a.size(); ++j)
1158 std::vector<std::int8_t> bc_markers1;
1159 std::vector<T> bc_values1;
1160 if (a[j] and !bcs1[j].empty())
1163 std::shared_ptr<const mesh::Mesh<U>> mesh = a[j]->get().mesh();
1165 throw std::runtime_error(
"Unable to extract a mesh.");
1166 mdspan2_t x_dofmap = mesh->geometry().dofmap();
1167 std::span _x = mesh->geometry().x();
1168 md::mdspan<const scalar_value_t<T>,
1169 md::extents<std::size_t, md::dynamic_extent, 3>>
1170 x(_x.data(), _x.size() / 3, 3);
1172 assert(a[j]->get().function_spaces().at(0));
1173 auto V1 = a[j]->get().function_spaces()[1];
1175 auto map1 = V1->dofmap()->index_map;
1176 const int bs1 = V1->dofmap()->index_map_bs();
1178 const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
1179 bc_markers1.assign(crange,
false);
1180 bc_values1.assign(crange, 0);
1181 for (
auto& bc : bcs1[j])
1183 bc.get().mark_dofs(bc_markers1);
1184 bc.get().set(bc_values1, std::nullopt, 1);
1189 lift_bc(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j],
1190 std::span<const T>(bc_values1), bc_markers1, x0[j], alpha);
1194 lift_bc(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j],
1195 std::span<const T>(bc_values1), bc_markers1,
1196 std::span<const T>(), alpha);
1209template <
typename V, std::floating_point U,
1210 dolfinx::scalar T =
typename std::remove_cvref_t<V>::value_type>
1211 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
1212void assemble_vector(
1213 V&& b,
const Form<T, U>& L,
1214 md::mdspan<
const scalar_value_t<T>,
1215 md::extents<std::size_t, md::dynamic_extent, 3>>
1217 std::span<const T> constants,
1218 const std::map<std::pair<IntegralType, int>,
1219 std::pair<std::span<const T>,
int>>& coefficients)
1222 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1226 auto mesh0 = L.function_spaces().at(0)->mesh();
1229 const int num_cell_types = mesh->topology()->cell_types().size();
1230 for (
int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
1233 mdspan2_t x_dofmap = mesh->geometry().dofmap(cell_type_idx);
1236 assert(L.function_spaces().at(0));
1237 auto element = L.function_spaces().at(0)->elements(cell_type_idx);
1239 std::shared_ptr<const fem::DofMap> dofmap
1240 = L.function_spaces().at(0)->dofmaps(cell_type_idx);
1242 auto dofs = dofmap->map();
1243 const int bs = dofmap->bs();
1245 fem::DofTransformKernel<T>
auto P0
1248 std::span<const std::uint32_t> cell_info0;
1249 if (element->needs_dof_transformations() or L.needs_facet_permutations())
1251 mesh0->topology_mutable()->create_entity_permutations();
1252 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
1262 assert(
cells.size() * cstride == coeffs.size());
1265 impl::assemble_cells<1>(
1266 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1267 md::mdspan(coeffs.data(),
cells.size(), cstride), cell_info0);
1271 impl::assemble_cells<3>(
1272 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1273 md::mdspan(coeffs.data(),
cells.size(), cstride), cell_info0);
1277 impl::assemble_cells(
1278 P0, b, x_dofmap, x, cells, {dofs, bs, cells0}, fn, constants,
1279 md::mdspan(coeffs.data(),
cells.size(), cstride), cell_info0);
1283 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
1284 if (L.needs_facet_permutations())
1286 mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx];
1287 int num_facets_per_cell
1289 mesh->topology_mutable()->create_entity_permutations();
1290 const std::vector<std::uint8_t>& p
1291 = mesh->topology()->get_facet_permutations();
1292 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
1293 num_facets_per_cell);
1297 = md::mdspan<
const std::int32_t,
1298 md::extents<std::size_t, md::dynamic_extent, 2>>;
1304 auto& [coeffs, cstride]
1307 mdspanx2_t facets(f.data(), f.size() / 2, 2);
1309 mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
1310 assert((facets.size() / 2) * cstride == coeffs.size());
1313 impl::assemble_exterior_facets<1>(
1314 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1315 md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0,
1320 impl::assemble_exterior_facets<3>(
1321 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1322 md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
1327 impl::assemble_exterior_facets(
1328 P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants,
1329 md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0,
1337 = md::mdspan<
const std::int32_t,
1338 md::extents<std::size_t, md::dynamic_extent, 2, 2>>;
1340 = md::mdspan<
const T, md::extents<std::size_t, md::dynamic_extent, 2,
1341 md::dynamic_extent>>;
1345 auto& [coeffs, cstride]
1349 assert((facets.size() / 4) * 2 * cstride == coeffs.size());
1352 impl::assemble_interior_facets<1>(
1354 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1356 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1358 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1363 impl::assemble_interior_facets<3>(
1365 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1367 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1369 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1374 impl::assemble_interior_facets(
1376 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
1378 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
1380 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride),
1393template <
typename V, std::floating_point U,
1394 dolfinx::scalar T =
typename std::remove_cvref_t<V>::value_type>
1395 requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
1396void assemble_vector(
1397 V&& b,
const Form<T, U>& L, std::span<const T> constants,
1398 const std::map<std::pair<IntegralType, int>,
1399 std::pair<std::span<const T>,
int>>& coefficients)
1402 = md::mdspan<const scalar_value_t<T>,
1403 md::extents<std::size_t, md::dynamic_extent, 3>>;
1405 std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
1407 auto x = mesh->geometry().x();
1408 if constexpr (std::is_same_v<U, scalar_value_t<T>>)
1410 impl::assemble_vector(b, L, mdspanx3_t(x.data(), x.size() / 3, 3),
1411 constants, coefficients);
1415 std::vector<scalar_value_t<T>> _x(x.begin(), x.end());
1416 impl::assemble_vector(b, L, mdspanx3_t(_x.data(), _x.size() / 3, 3),
1417 constants, coefficients);
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: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