DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
assemble_scalar_impl.h
1// Copyright (C) 2019-2025 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 "Form.h"
11#include "FunctionSpace.h"
12#include "utils.h"
13#include <algorithm>
14#include <basix/mdspan.hpp>
15#include <dolfinx/common/IndexMap.h>
16#include <dolfinx/mesh/Geometry.h>
17#include <dolfinx/mesh/Mesh.h>
18#include <dolfinx/mesh/Topology.h>
19#include <memory>
20#include <vector>
21
22namespace dolfinx::fem::impl
23{
25template <dolfinx::scalar T>
26T assemble_cells(mdspan2_t x_dofmap,
27 md::mdspan<const scalar_value_t<T>,
28 md::extents<std::size_t, md::dynamic_extent, 3>>
29 x,
30 std::span<const std::int32_t> cells, FEkernel<T> auto fn,
31 std::span<const T> constants,
32 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs)
33{
34 T value(0);
35 if (cells.empty())
36 return value;
37
38 // Create data structures used in assembly
39 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
40
41 // Iterate over all cells
42 for (std::size_t index = 0; index < cells.size(); ++index)
43 {
44 std::int32_t c = cells[index];
45
46 // Get cell coordinates/geometry
47 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
48 for (std::size_t i = 0; i < x_dofs.size(); ++i)
49 {
50 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
51 }
52
53 fn(&value, &coeffs(index, 0), constants.data(), cdofs.data(), nullptr,
54 nullptr, nullptr);
55 }
56
57 return value;
58}
59
61template <dolfinx::scalar T>
62T assemble_exterior_facets(
63 mdspan2_t x_dofmap,
64 md::mdspan<const scalar_value_t<T>,
65 md::extents<std::size_t, md::dynamic_extent, 3>>
66 x,
67 md::mdspan<const std::int32_t,
68 md::extents<std::size_t, md::dynamic_extent, 2>>
69 facets,
70 FEkernel<T> auto fn, std::span<const T> constants,
71 md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
72 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
73{
74 T value(0);
75 if (facets.empty())
76 return value;
77
78 // Create data structures used in assembly
79 std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
80
81 // Iterate over all facets
82 for (std::size_t f = 0; f < facets.extent(0); ++f)
83 {
84 std::int32_t cell = facets(f, 0);
85 std::int32_t local_facet = facets(f, 1);
86
87 // Get cell coordinates/geometry
88 auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
89 for (std::size_t i = 0; i < x_dofs.size(); ++i)
90 std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
91
92 // Permutations
93 std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet);
94 fn(&value, &coeffs(f, 0), constants.data(), cdofs.data(), &local_facet,
95 &perm, nullptr);
96 }
97
98 return value;
99}
100
102template <dolfinx::scalar T>
103T assemble_interior_facets(
104 mdspan2_t x_dofmap,
105 md::mdspan<const scalar_value_t<T>,
106 md::extents<std::size_t, md::dynamic_extent, 3>>
107 x,
108 md::mdspan<const std::int32_t,
109 md::extents<std::size_t, md::dynamic_extent, 2, 2>>
110 facets,
111 FEkernel<T> auto fn, std::span<const T> constants,
112 md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
113 md::dynamic_extent>>
114 coeffs,
115 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms)
116{
117 T value(0);
118 if (facets.empty())
119 return value;
120
121 // Create data structures used in assembly
122 using X = scalar_value_t<T>;
123 std::vector<X> cdofs(2 * x_dofmap.extent(1) * 3);
124 std::span<X> cdofs0(cdofs.data(), x_dofmap.extent(1) * 3);
125 std::span<X> cdofs1(cdofs.data() + x_dofmap.extent(1) * 3,
126 x_dofmap.extent(1) * 3);
127
128 // Iterate over all facets
129 for (std::size_t f = 0; f < facets.extent(0); ++f)
130 {
131 std::array cells = {facets(f, 0, 0), facets(f, 1, 0)};
132 std::array local_facet = {facets(f, 0, 1), facets(f, 1, 1)};
133
134 // Get cell geometry
135 auto x_dofs0 = md::submdspan(x_dofmap, cells[0], md::full_extent);
136 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
137 std::copy_n(&x(x_dofs0[i], 0), 3, std::next(cdofs0.begin(), 3 * i));
138 auto x_dofs1 = md::submdspan(x_dofmap, cells[1], md::full_extent);
139 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
140 std::copy_n(&x(x_dofs1[i], 0), 3, std::next(cdofs1.begin(), 3 * i));
141
142 std::array perm = perms.empty()
143 ? std::array<std::uint8_t, 2>{0, 0}
144 : std::array{perms(cells[0], local_facet[0]),
145 perms(cells[1], local_facet[1])};
146 fn(&value, &coeffs(f, 0, 0), constants.data(), cdofs.data(),
147 local_facet.data(), perm.data(), nullptr);
148 }
149
150 return value;
151}
152
154template <dolfinx::scalar T, std::floating_point U>
155T assemble_scalar(
156 const fem::Form<T, U>& M, mdspan2_t x_dofmap,
157 md::mdspan<const scalar_value_t<T>,
158 md::extents<std::size_t, md::dynamic_extent, 3>>
159 x,
160 std::span<const T> constants,
161 const std::map<std::pair<IntegralType, int>,
162 std::pair<std::span<const T>, int>>& coefficients)
163{
164 std::shared_ptr<const mesh::Mesh<U>> mesh = M.mesh();
165 assert(mesh);
166
167 T value = 0;
168 for (int i : M.integral_ids(IntegralType::cell))
169 {
170 auto fn = M.kernel(IntegralType::cell, i, 0);
171 assert(fn);
172 auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
173 std::span<const std::int32_t> cells = M.domain(IntegralType::cell, i, 0);
174 assert(cells.size() * cstride == coeffs.size());
175 value += impl::assemble_cells(
176 x_dofmap, x, cells, fn, constants,
177 md::mdspan(coeffs.data(), cells.size(), cstride));
178 }
179
180 mesh::CellType cell_type = mesh->topology()->cell_type();
181 int num_facets_per_cell
182 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
183 md::mdspan<const std::uint8_t, md::dextents<std::size_t, 2>> perms;
184 if (M.needs_facet_permutations())
185 {
186 mesh->topology_mutable()->create_entity_permutations();
187 const std::vector<std::uint8_t>& p
188 = mesh->topology()->get_facet_permutations();
189 perms = md::mdspan(p.data(), p.size() / num_facets_per_cell,
190 num_facets_per_cell);
191 }
192
193 for (int i : M.integral_ids(IntegralType::exterior_facet))
194 {
195 auto fn = M.kernel(IntegralType::exterior_facet, i, 0);
196 assert(fn);
197 auto& [coeffs, cstride]
198 = coefficients.at({IntegralType::exterior_facet, i});
199
200 std::span facets = M.domain(IntegralType::exterior_facet, i, 0);
201 assert((facets.size() / 2) * cstride == coeffs.size());
202 value += impl::assemble_exterior_facets(
203 x_dofmap, x,
204 md::mdspan<const std::int32_t,
205 md::extents<std::size_t, md::dynamic_extent, 2>>(
206 facets.data(), facets.size() / 2, 2),
207 fn, constants, md::mdspan(coeffs.data(), facets.size() / 2, cstride),
208 perms);
209 }
210
211 for (int i : M.integral_ids(IntegralType::interior_facet))
212 {
213 auto fn = M.kernel(IntegralType::interior_facet, i, 0);
214 assert(fn);
215 auto& [coeffs, cstride]
216 = coefficients.at({IntegralType::interior_facet, i});
217 std::span facets = M.domain(IntegralType::interior_facet, i, 0);
218 assert((facets.size() / 4) * 2 * cstride == coeffs.size());
219 value += impl::assemble_interior_facets(
220 x_dofmap, x,
221 md::mdspan<const std::int32_t,
222 md::extents<std::size_t, md::dynamic_extent, 2, 2>>(
223 facets.data(), facets.size() / 4, 2, 2),
224 fn, constants,
225 md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
226 md::dynamic_extent>>(
227 coeffs.data(), facets.size() / 4, 2, cstride),
228 perms);
229 }
230
231 return value;
232}
233
234} // namespace dolfinx::fem::impl
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.
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