11 #include <dolfinx/common/IndexMap.h>
12 #include <dolfinx/common/types.h>
13 #include <dolfinx/fem/Constant.h>
14 #include <dolfinx/fem/FunctionSpace.h>
15 #include <dolfinx/mesh/Geometry.h>
16 #include <dolfinx/mesh/Mesh.h>
17 #include <dolfinx/mesh/Topology.h>
21 namespace dolfinx::fem::impl
31 const mesh::Geometry& geometry,
32 const std::vector<std::int32_t>& active_cells,
33 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
34 const std::uint8_t*,
const std::uint32_t)>& fn,
35 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
36 const std::vector<std::uint32_t>& cell_info);
40 T assemble_exterior_facets(
41 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
42 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
43 const std::uint8_t*,
const std::uint32_t)>& fn,
44 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
45 const std::vector<std::uint32_t>& cell_info,
46 const std::vector<std::uint8_t>& perms);
50 T assemble_interior_facets(
51 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
52 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
53 const std::uint8_t*,
const std::uint32_t)>& fn,
54 const array2d<T>& coeffs,
const std::vector<int>& offsets,
55 const std::vector<T>& constant_values,
56 const std::vector<std::uint32_t>& cell_info,
57 const std::vector<std::uint8_t>& perms);
63 std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
65 const int tdim = mesh->topology().dim();
66 const std::int32_t num_cells
67 = mesh->topology().connectivity(tdim, 0)->num_nodes();
70 const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants
73 std::vector<T> constant_values;
74 for (
auto const& constant : constants)
77 const std::vector<T>& array = constant->value;
78 constant_values.insert(constant_values.end(), array.begin(), array.end());
84 const bool needs_permutation_data = M.needs_permutation_data();
85 if (needs_permutation_data)
86 mesh->topology_mutable().create_entity_permutations();
87 const std::vector<std::uint32_t>& cell_info
88 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
89 : std::vector<std::uint32_t>(num_cells);
92 for (
int i : M.integral_ids(IntegralType::cell))
94 const auto& fn = M.kernel(IntegralType::cell, i);
95 const std::vector<std::int32_t>& active_cells
96 = M.domains(IntegralType::cell, i);
97 value += fem::impl::assemble_cells(mesh->geometry(), active_cells, fn,
98 coeffs, constant_values, cell_info);
101 if (M.num_integrals(IntegralType::exterior_facet) > 0
102 or M.num_integrals(IntegralType::interior_facet) > 0)
105 mesh->topology_mutable().create_entities(tdim - 1);
106 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
107 mesh->topology_mutable().create_entity_permutations();
109 const std::vector<std::uint8_t>& perms
110 = mesh->topology().get_facet_permutations();
112 for (
int i : M.integral_ids(IntegralType::exterior_facet))
114 const auto& fn = M.kernel(IntegralType::exterior_facet, i);
115 const std::vector<std::int32_t>& active_facets
116 = M.domains(IntegralType::exterior_facet, i);
117 value += fem::impl::assemble_exterior_facets(
118 *mesh, active_facets, fn, coeffs, constant_values, cell_info, perms);
121 const std::vector<int> c_offsets = M.coefficient_offsets();
122 for (
int i : M.integral_ids(IntegralType::interior_facet))
124 const auto& fn = M.kernel(IntegralType::interior_facet, i);
125 const std::vector<std::int32_t>& active_facets
126 = M.domains(IntegralType::interior_facet, i);
127 value += fem::impl::assemble_interior_facets(
128 *mesh, active_facets, fn, coeffs, c_offsets, constant_values,
136 template <
typename T>
138 const mesh::Geometry& geometry,
139 const std::vector<std::int32_t>& active_cells,
140 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
141 const std::uint8_t*,
const std::uint32_t)>& fn,
142 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
143 const std::vector<std::uint32_t>& cell_info)
145 const int gdim = geometry.dim();
148 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
151 const int num_dofs_g = x_dofmap.num_links(0);
152 const array2d<double>& x_g = geometry.x();
155 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
159 for (std::int32_t c : active_cells)
162 auto x_dofs = x_dofmap.links(c);
163 for (std::size_t i = 0; i < x_dofs.size(); ++i)
165 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
166 std::next(coordinate_dofs.begin(), i * gdim));
169 auto coeff_cell = coeffs.row(c);
170 fn(&value, coeff_cell.data(), constant_values.data(),
171 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
177 template <
typename T>
178 T assemble_exterior_facets(
179 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
180 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
181 const std::uint8_t*,
const std::uint32_t)>& fn,
182 const array2d<T>& coeffs,
const std::vector<T>& constant_values,
183 const std::vector<std::uint32_t>& cell_info,
184 const std::vector<std::uint8_t>& perms)
186 const int gdim = mesh.geometry().dim();
187 const int tdim = mesh.topology().dim();
190 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
193 const int num_dofs_g = x_dofmap.num_links(0);
194 const array2d<double>& x_g = mesh.geometry().x();
197 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
199 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
201 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
206 for (std::int32_t facet : active_facets)
209 assert(f_to_c->num_links(facet) == 1);
210 const int cell = f_to_c->links(facet)[0];
213 auto facets = c_to_f->links(cell);
214 auto it = std::find(facets.begin(), facets.end(), facet);
215 assert(it != facets.end());
216 const int local_facet = std::distance(facets.data(), it);
219 auto x_dofs = x_dofmap.links(cell);
220 for (std::size_t i = 0; i < x_dofs.size(); ++i)
222 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
223 std::next(coordinate_dofs.begin(), i * gdim));
226 auto coeff_cell = coeffs.row(cell);
227 fn(&value, coeff_cell.data(), constant_values.data(),
228 coordinate_dofs.data(), &local_facet,
229 &perms[cell * facets.size() + local_facet], cell_info[cell]);
235 template <
typename T>
236 T assemble_interior_facets(
237 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
238 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
239 const std::uint8_t*,
const std::uint32_t)>& fn,
240 const array2d<T>& coeffs,
const std::vector<int>& offsets,
241 const std::vector<T>& constant_values,
242 const std::vector<std::uint32_t>& cell_info,
243 const std::vector<std::uint8_t>& perms)
245 const int gdim = mesh.geometry().dim();
246 const int tdim = mesh.topology().dim();
249 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
252 const int num_dofs_g = x_dofmap.num_links(0);
253 const array2d<double>& x_g = mesh.geometry().x();
256 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
257 std::vector<T> coeff_array(2 * offsets.back());
258 assert(offsets.back() == coeffs.shape[1]);
260 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
262 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
267 const int offset_g = gdim * num_dofs_g;
268 for (std::int32_t f : active_facets)
271 auto cells = f_to_c->links(f);
272 assert(
cells.size() == 2);
274 const int facets_per_cell = c_to_f->num_links(cells[0]);
277 std::array<int, 2> local_facet;
278 for (
int i = 0; i < 2; ++i)
280 auto facets = c_to_f->links(cells[i]);
281 auto it = std::find(facets.begin(), facets.end(), f);
282 assert(it != facets.end());
283 local_facet[i] = std::distance(facets.begin(), it);
287 auto x_dofs0 = x_dofmap.links(cells[0]);
288 auto x_dofs1 = x_dofmap.links(cells[1]);
289 for (
int i = 0; i < num_dofs_g; ++i)
291 for (
int j = 0; j < gdim; ++j)
293 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
294 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
300 auto coeff_cell0 = coeffs.row(cells[0]);
301 auto coeff_cell1 = coeffs.row(cells[1]);
304 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
307 const int num_entries = offsets[i + 1] - offsets[i];
308 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
309 std::next(coeff_array.begin(), 2 * offsets[i]));
310 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
311 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
314 const std::array perm{perms[
cells[0] * facets_per_cell + local_facet[0]],
315 perms[
cells[1] * facets_per_cell + local_facet[1]]};
316 fn(&value, coeff_array.data(), constant_values.data(),
317 coordinate_dofs.data(), local_facet.data(), perm.data(),
318 cell_info[cells[0]]);
void cells(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const std::reference_wrapper< const fem::DofMap >, 2 > &dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:18
T assemble_scalar(const Form< T > &M)
Assemble functional into scalar. Caller is responsible for accumulation across processes.
Definition: assembler.h:33
array2d< typename U::scalar_type > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:398