11#include "FunctionSpace.h"
15#include <dolfinx/la/utils.h>
16#include <dolfinx/mesh/Geometry.h>
17#include <dolfinx/mesh/Mesh.h>
18#include <dolfinx/mesh/Topology.h>
25namespace dolfinx::fem::impl
28using mdspan2_t = md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>>;
63template <dolfinx::scalar T>
65 la::MatSet<T>
auto mat_set, mdspan2_t x_dofmap,
66 md::mdspan<
const scalar_value_t<T>,
67 md::extents<std::size_t, md::dynamic_extent, 3>>
69 std::span<const std::int32_t> cells,
70 std::tuple<mdspan2_t,
int, std::span<const std::int32_t>> dofmap0,
71 fem::DofTransformKernel<T>
auto P0,
72 std::tuple<mdspan2_t,
int, std::span<const std::int32_t>> dofmap1,
73 fem::DofTransformKernel<T>
auto P1T, std::span<const std::int8_t> bc0,
74 std::span<const std::int8_t> bc1, FEkernel<T>
auto kernel,
75 md::mdspan<
const T, md::dextents<std::size_t, 2>> coeffs,
76 std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
77 std::span<const std::uint32_t> cell_info1)
82 const auto [dmap0, bs0, cells0] = dofmap0;
83 const auto [dmap1, bs1, cells1] = dofmap1;
86 const int num_dofs0 = dmap0.extent(1);
87 const int num_dofs1 = dmap1.extent(1);
88 const int ndim0 = bs0 * num_dofs0;
89 const int ndim1 = bs1 * num_dofs1;
90 std::vector<T> Ae(ndim0 * ndim1);
92 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
95 assert(cells0.size() ==
cells.size());
96 assert(cells1.size() ==
cells.size());
97 for (std::size_t c = 0; c <
cells.size(); ++c)
102 std::int32_t cell0 = cells0[c];
103 std::int32_t cell1 = cells1[c];
106 auto x_dofs = md::submdspan(x_dofmap,
cell, md::full_extent);
107 for (std::size_t i = 0; i < x_dofs.size(); ++i)
108 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
111 std::ranges::fill(Ae, 0);
112 kernel(Ae.data(), &coeffs(c, 0), constants.data(), cdofs.data(),
nullptr,
116 P0(_Ae, cell_info0, cell0, ndim1);
117 P1T(_Ae, cell_info1, cell1, ndim0);
120 std::span dofs0(dmap0.data_handle() + cell0 * num_dofs0, num_dofs0);
121 std::span dofs1(dmap1.data_handle() + cell1 * num_dofs1, num_dofs1);
125 for (
int i = 0; i < num_dofs0; ++i)
127 for (
int k = 0; k < bs0; ++k)
129 if (bc0[bs0 * dofs0[i] + k])
132 const int row = bs0 * i + k;
133 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0);
141 for (
int j = 0; j < num_dofs1; ++j)
143 for (
int k = 0; k < bs1; ++k)
145 if (bc1[bs1 * dofs1[j] + k])
148 const int col = bs1 * j + k;
149 for (
int row = 0; row < ndim0; ++row)
150 Ae[row * ndim1 + col] = 0;
156 mat_set(dofs0, dofs1, Ae);
194template <dolfinx::scalar T>
195void assemble_exterior_facets(
196 la::MatSet<T>
auto mat_set, mdspan2_t x_dofmap,
197 md::mdspan<
const scalar_value_t<T>,
198 md::extents<std::size_t, md::dynamic_extent, 3>>
200 md::mdspan<
const std::int32_t,
201 std::extents<std::size_t, md::dynamic_extent, 2>>
203 std::tuple<mdspan2_t,
int,
204 md::mdspan<
const std::int32_t,
205 std::extents<std::size_t, md::dynamic_extent, 2>>>
207 fem::DofTransformKernel<T>
auto P0,
208 std::tuple<mdspan2_t,
int,
209 md::mdspan<
const std::int32_t,
210 std::extents<std::size_t, md::dynamic_extent, 2>>>
212 fem::DofTransformKernel<T>
auto P1T, std::span<const std::int8_t> bc0,
213 std::span<const std::int8_t> bc1, FEkernel<T>
auto kernel,
214 md::mdspan<
const T, md::dextents<std::size_t, 2>> coeffs,
215 std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
216 std::span<const std::uint32_t> cell_info1,
217 md::mdspan<
const std::uint8_t, md::dextents<std::size_t, 2>> perms)
222 const auto [dmap0, bs0, facets0] = dofmap0;
223 const auto [dmap1, bs1, facets1] = dofmap1;
226 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
227 const int num_dofs0 = dmap0.extent(1);
228 const int num_dofs1 = dmap1.extent(1);
229 const int ndim0 = bs0 * num_dofs0;
230 const int ndim1 = bs1 * num_dofs1;
231 std::vector<T> Ae(ndim0 * ndim1);
232 std::span<T> _Ae(Ae);
233 assert(facets0.size() == facets.size());
234 assert(facets1.size() == facets.size());
235 for (std::size_t f = 0; f < facets.extent(0); ++f)
240 std::int32_t
cell = facets(f, 0);
241 std::int32_t local_facet = facets(f, 1);
242 std::int32_t cell0 = facets0(f, 0);
243 std::int32_t cell1 = facets1(f, 0);
246 auto x_dofs = md::submdspan(x_dofmap,
cell, md::full_extent);
247 for (std::size_t i = 0; i < x_dofs.size(); ++i)
248 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
251 std::uint8_t perm = perms.empty() ? 0 : perms(
cell, local_facet);
254 std::ranges::fill(Ae, 0);
255 kernel(Ae.data(), &coeffs(f, 0), constants.data(), cdofs.data(),
256 &local_facet, &perm,
nullptr);
258 P0(_Ae, cell_info0, cell0, ndim1);
259 P1T(_Ae, cell_info1, cell1, ndim0);
262 std::span dofs0(dmap0.data_handle() + cell0 * num_dofs0, num_dofs0);
263 std::span dofs1(dmap1.data_handle() + cell1 * num_dofs1, num_dofs1);
266 for (
int i = 0; i < num_dofs0; ++i)
268 for (
int k = 0; k < bs0; ++k)
270 if (bc0[bs0 * dofs0[i] + k])
273 const int row = bs0 * i + k;
274 std::fill_n(std::next(Ae.begin(), ndim1 * row), ndim1, 0);
281 for (
int j = 0; j < num_dofs1; ++j)
283 for (
int k = 0; k < bs1; ++k)
285 if (bc1[bs1 * dofs1[j] + k])
288 const int col = bs1 * j + k;
289 for (
int row = 0; row < ndim0; ++row)
290 Ae[row * ndim1 + col] = 0;
296 mat_set(dofs0, dofs1, Ae);
334template <dolfinx::scalar T>
335void assemble_interior_facets(
336 la::MatSet<T>
auto mat_set, mdspan2_t x_dofmap,
337 md::mdspan<
const scalar_value_t<T>,
338 md::extents<std::size_t, md::dynamic_extent, 3>>
340 md::mdspan<
const std::int32_t,
341 std::extents<std::size_t, md::dynamic_extent, 2, 2>>
343 std::tuple<
const DofMap&,
int,
344 md::mdspan<
const std::int32_t,
345 std::extents<std::size_t, md::dynamic_extent, 2, 2>>>
347 fem::DofTransformKernel<T>
auto P0,
348 std::tuple<
const DofMap&,
int,
349 md::mdspan<
const std::int32_t,
350 std::extents<std::size_t, md::dynamic_extent, 2, 2>>>
352 fem::DofTransformKernel<T>
auto P1T, std::span<const std::int8_t> bc0,
353 std::span<const std::int8_t> bc1, FEkernel<T>
auto kernel,
354 md::mdspan<
const T, md::extents<std::size_t, md::dynamic_extent, 2,
357 std::span<const T> constants, std::span<const std::uint32_t> cell_info0,
358 std::span<const std::uint32_t> cell_info1,
359 md::mdspan<
const std::uint8_t, md::dextents<std::size_t, 2>> perms)
364 const auto [dmap0, bs0, facets0] = dofmap0;
365 const auto [dmap1, bs1, facets1] = dofmap1;
368 using X = scalar_value_t<T>;
369 std::vector<X> cdofs(2 * x_dofmap.extent(1) * 3);
370 std::span<X> cdofs0(cdofs.data(), x_dofmap.extent(1) * 3);
371 std::span<X> cdofs1(cdofs.data() + x_dofmap.extent(1) * 3,
372 x_dofmap.extent(1) * 3);
374 std::vector<T> Ae, be;
377 std::vector<std::int32_t> dmapjoint0, dmapjoint1;
378 assert(facets0.size() == facets.size());
379 assert(facets1.size() == facets.size());
380 for (std::size_t f = 0; f < facets.extent(0); ++f)
384 std::array
cells{facets(f, 0, 0), facets(f, 1, 0)};
385 std::array cells0{facets0(f, 0, 0), facets0(f, 1, 0)};
386 std::array cells1{facets1(f, 0, 0), facets1(f, 1, 0)};
389 std::array local_facet{facets(f, 0, 1), facets(f, 1, 1)};
392 auto x_dofs0 = md::submdspan(x_dofmap, cells[0], md::full_extent);
393 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
394 std::copy_n(&x(x_dofs0[i], 0), 3, std::next(cdofs0.begin(), 3 * i));
395 auto x_dofs1 = md::submdspan(x_dofmap, cells[1], md::full_extent);
396 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
397 std::copy_n(&x(x_dofs1[i], 0), 3, std::next(cdofs1.begin(), 3 * i));
400 std::span<const std::int32_t> dmap0_cell0 = dmap0.cell_dofs(cells0[0]);
401 std::span<const std::int32_t> dmap0_cell1 = dmap0.cell_dofs(cells0[1]);
402 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
403 std::ranges::copy(dmap0_cell0, dmapjoint0.begin());
404 std::ranges::copy(dmap0_cell1,
405 std::next(dmapjoint0.begin(), dmap0_cell0.size()));
407 std::span<const std::int32_t> dmap1_cell0 = dmap1.cell_dofs(cells1[0]);
408 std::span<const std::int32_t> dmap1_cell1 = dmap1.cell_dofs(cells1[1]);
409 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
410 std::ranges::copy(dmap1_cell0, dmapjoint1.begin());
411 std::ranges::copy(dmap1_cell1,
412 std::next(dmapjoint1.begin(), dmap1_cell0.size()));
414 const int num_rows = bs0 * dmapjoint0.size();
415 const int num_cols = bs1 * dmapjoint1.size();
418 Ae.resize(num_rows * num_cols);
419 std::ranges::fill(Ae, 0);
421 std::array perm = perms.empty()
422 ? std::array<std::uint8_t, 2>{0, 0}
423 : std::array{perms(cells[0], local_facet[0]),
424 perms(cells[1], local_facet[1])};
425 kernel(Ae.data(), &coeffs(f, 0, 0), constants.data(), cdofs.data(),
426 local_facet.data(), perm.data(),
nullptr);
435 std::span<T> _Ae(Ae);
436 std::span<T> sub_Ae0 = _Ae.subspan(bs0 * dmap0_cell0.size() * num_cols,
437 bs0 * dmap0_cell1.size() * num_cols);
439 P0(_Ae, cell_info0, cells0[0], num_cols);
440 P0(sub_Ae0, cell_info0, cells0[1], num_cols);
441 P1T(_Ae, cell_info1, cells1[0], num_rows);
443 for (
int row = 0; row < num_rows; ++row)
447 std::span<T> sub_Ae1 = _Ae.subspan(
448 row * num_cols + bs1 * dmap1_cell0.size(), bs1 * dmap1_cell1.size());
449 P1T(sub_Ae1, cell_info1, cells1[1], 1);
455 for (std::size_t i = 0; i < dmapjoint0.size(); ++i)
457 for (
int k = 0; k < bs0; ++k)
459 if (bc0[bs0 * dmapjoint0[i] + k])
462 std::fill_n(std::next(Ae.begin(), num_cols * (bs0 * i + k)),
470 for (std::size_t j = 0; j < dmapjoint1.size(); ++j)
472 for (
int k = 0; k < bs1; ++k)
474 if (bc1[bs1 * dmapjoint1[j] + k])
477 for (
int m = 0; m < num_rows; ++m)
478 Ae[m * num_cols + bs1 * j + k] = 0;
484 mat_set(dmapjoint0, dmapjoint1, Ae);
512template <dolfinx::scalar T, std::
floating_po
int U>
514 la::MatSet<T>
auto mat_set,
const Form<T, U>& a,
515 md::mdspan<
const scalar_value_t<T>,
516 md::extents<std::size_t, md::dynamic_extent, 3>>
518 std::span<const T> constants,
519 const std::map<std::pair<IntegralType, int>,
520 std::pair<std::span<const T>,
int>>& coefficients,
521 std::span<const std::int8_t> bc0, std::span<const std::int8_t> bc1)
524 std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
528 auto mesh0 = a.function_spaces().at(0)->mesh();
532 auto mesh1 = a.function_spaces().at(1)->mesh();
541 const int num_cell_types = mesh->topology()->cell_types().size();
542 for (
int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
545 mdspan2_t x_dofmap = mesh->geometry().dofmap(cell_type_idx);
548 std::shared_ptr<const fem::DofMap> dofmap0
549 = a.function_spaces().at(0)->dofmaps(cell_type_idx);
550 std::shared_ptr<const fem::DofMap> dofmap1
551 = a.function_spaces().at(1)->dofmaps(cell_type_idx);
554 auto dofs0 = dofmap0->map();
555 const int bs0 = dofmap0->bs();
556 auto dofs1 = dofmap1->map();
557 const int bs1 = dofmap1->bs();
559 auto element0 = a.function_spaces().at(0)->elements(cell_type_idx);
561 auto element1 = a.function_spaces().at(1)->elements(cell_type_idx);
563 fem::DofTransformKernel<T>
auto P0
565 fem::DofTransformKernel<T>
auto P1T
566 = element1->template dof_transformation_right_fn<T>(
569 std::span<const std::uint32_t> cell_info0;
570 std::span<const std::uint32_t> cell_info1;
571 if (element0->needs_dof_transformations()
572 or element1->needs_dof_transformations()
573 or a.needs_facet_permutations())
575 mesh0->topology_mutable()->create_entity_permutations();
576 mesh1->topology_mutable()->create_entity_permutations();
577 cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info());
578 cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info());
589 assert(
cells.size() * cstride == coeffs.size());
590 impl::assemble_cells(mat_set, x_dofmap, x, cells, {dofs0, bs0, cells0},
591 P0, {dofs1, bs1, cells1}, P1T, bc0, bc1, fn,
592 md::mdspan(coeffs.data(),
cells.size(), cstride),
593 constants, cell_info0, cell_info1);
596 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
597 if (a.needs_facet_permutations())
599 mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx];
600 int num_facets_per_cell
602 mesh->topology_mutable()->create_entity_permutations();
603 const std::vector<std::uint8_t>& p
604 = mesh->topology()->get_facet_permutations();
605 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
606 num_facets_per_cell);
611 if (num_cell_types > 1)
613 throw std::runtime_error(
"Exterior facet integrals with mixed "
614 "topology aren't supported yet");
618 = md::mdspan<
const std::int32_t,
619 md::extents<std::size_t, md::dynamic_extent, 2>>;
623 auto& [coeffs, cstride]
627 mdspanx2_t facets(f.data(), f.size() / 2, 2);
629 mdspanx2_t facets0(f0.data(), f0.size() / 2, 2);
631 mdspanx2_t facets1(f1.data(), f1.size() / 2, 2);
632 assert((facets.size() / 2) * cstride == coeffs.size());
633 impl::assemble_exterior_facets(
634 mat_set, x_dofmap, x, facets, {dofs0, bs0, facets0}, P0,
635 {dofs1, bs1, facets1}, P1T, bc0, bc1, fn,
636 md::mdspan(coeffs.data(), facets.extent(0), cstride), constants,
637 cell_info0, cell_info1, perms);
642 if (num_cell_types > 1)
644 throw std::runtime_error(
"Interior facet integrals with mixed "
645 "topology aren't supported yet");
649 = md::mdspan<
const std::int32_t,
650 md::extents<std::size_t, md::dynamic_extent, 2, 2>>;
652 = md::mdspan<
const T, md::extents<std::size_t, md::dynamic_extent, 2,
653 md::dynamic_extent>>;
657 auto& [coeffs, cstride]
663 assert((facets.size() / 4) * 2 * cstride == coeffs.size());
664 impl::assemble_interior_facets(
665 mat_set, x_dofmap, x,
666 mdspanx22_t(facets.data(), facets.size() / 4, 2, 2),
668 mdspanx22_t(facets0.data(), facets0.size() / 4, 2, 2)},
671 mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)},
673 mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), constants,
674 cell_info0, cell_info1, perms);
Degree-of-freedom map representations and tools.
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
@ 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