12 #include <dolfinx/common/IndexMap.h>
13 #include <dolfinx/common/types.h>
14 #include <dolfinx/fem/Constant.h>
15 #include <dolfinx/fem/FunctionSpace.h>
16 #include <dolfinx/mesh/Geometry.h>
17 #include <dolfinx/mesh/Mesh.h>
18 #include <dolfinx/mesh/Topology.h>
22 namespace dolfinx::fem::impl
32 const mesh::Geometry& geometry,
33 const std::vector<std::int32_t>& active_cells,
34 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
35 const std::uint8_t*,
const std::uint32_t)>& fn,
36 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
38 const std::vector<T>& constant_values,
39 const std::vector<std::uint32_t>& cell_info);
43 T assemble_exterior_facets(
44 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
45 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
46 const std::uint8_t*,
const std::uint32_t)>& fn,
47 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
49 const std::vector<T>& constant_values,
50 const std::vector<std::uint32_t>& cell_info,
51 const std::vector<std::uint8_t>& perms);
55 T assemble_interior_facets(
56 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
57 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
58 const std::uint8_t*,
const std::uint32_t)>& fn,
59 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
61 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
62 const std::vector<std::uint32_t>& cell_info,
63 const std::vector<std::uint8_t>& perms);
69 std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
71 const int tdim = mesh->topology().dim();
72 const std::int32_t num_cells
73 = mesh->topology().connectivity(tdim, 0)->num_nodes();
76 const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants
79 std::vector<T> constant_values;
80 for (
auto const& constant : constants)
83 const std::vector<T>& array = constant->value;
84 constant_values.insert(constant_values.end(), array.begin(), array.end());
88 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
91 const bool needs_permutation_data = M.needs_permutation_data();
92 if (needs_permutation_data)
93 mesh->topology_mutable().create_entity_permutations();
94 const std::vector<std::uint32_t>& cell_info
95 = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
96 : std::vector<std::uint32_t>(num_cells);
99 for (
int i : M.integral_ids(IntegralType::cell))
101 const auto& fn = M.kernel(IntegralType::cell, i);
102 const std::vector<std::int32_t>& active_cells
103 = M.domains(IntegralType::cell, i);
104 value += fem::impl::assemble_cells(mesh->geometry(), active_cells, fn,
105 coeffs, constant_values, cell_info);
108 if (M.num_integrals(IntegralType::exterior_facet) > 0
109 or M.num_integrals(IntegralType::interior_facet) > 0)
112 mesh->topology_mutable().create_entities(tdim - 1);
113 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
114 mesh->topology_mutable().create_entity_permutations();
116 const std::vector<std::uint8_t>& perms
117 = mesh->topology().get_facet_permutations();
119 for (
int i : M.integral_ids(IntegralType::exterior_facet))
121 const auto& fn = M.kernel(IntegralType::exterior_facet, i);
122 const std::vector<std::int32_t>& active_facets
123 = M.domains(IntegralType::exterior_facet, i);
124 value += fem::impl::assemble_exterior_facets(
125 *mesh, active_facets, fn, coeffs, constant_values, cell_info, perms);
128 const std::vector<int> c_offsets = M.coefficient_offsets();
129 for (
int i : M.integral_ids(IntegralType::interior_facet))
131 const auto& fn = M.kernel(IntegralType::interior_facet, i);
132 const std::vector<std::int32_t>& active_facets
133 = M.domains(IntegralType::interior_facet, i);
134 value += fem::impl::assemble_interior_facets(
135 *mesh, active_facets, fn, coeffs, c_offsets, constant_values,
143 template <
typename T>
145 const mesh::Geometry& geometry,
146 const std::vector<std::int32_t>& active_cells,
147 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
148 const std::uint8_t*,
const std::uint32_t)>& fn,
149 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
151 const std::vector<T>& constant_values,
152 const std::vector<std::uint32_t>& cell_info)
154 const int gdim = geometry.dim();
157 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
160 const int num_dofs_g = x_dofmap.num_links(0);
161 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
165 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
169 for (std::int32_t c : active_cells)
172 auto x_dofs = x_dofmap.links(c);
173 for (std::size_t i = 0; i < x_dofs.size(); ++i)
175 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
176 std::next(coordinate_dofs.begin(), i * gdim));
179 auto coeff_cell = coeffs.row(c);
180 fn(&value, coeff_cell.data(), constant_values.data(),
181 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
187 template <
typename T>
188 T assemble_exterior_facets(
189 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
190 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
191 const std::uint8_t*,
const std::uint32_t)>& fn,
192 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
194 const std::vector<T>& constant_values,
195 const std::vector<std::uint32_t>& cell_info,
196 const std::vector<std::uint8_t>& perms)
198 const int gdim = mesh.geometry().dim();
199 const int tdim = mesh.topology().dim();
202 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
205 const int num_dofs_g = x_dofmap.num_links(0);
206 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
207 = mesh.geometry().x();
210 std::vector<double> coordinate_dofs(num_dofs_g * gdim);
212 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
214 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
219 for (std::int32_t facet : active_facets)
222 assert(f_to_c->num_links(facet) == 1);
223 const int cell = f_to_c->links(facet)[0];
226 auto facets = c_to_f->links(cell);
227 auto it = std::find(facets.begin(), facets.end(), facet);
228 assert(it != facets.end());
229 const int local_facet = std::distance(facets.data(), it);
232 auto x_dofs = x_dofmap.links(cell);
233 for (std::size_t i = 0; i < x_dofs.size(); ++i)
235 std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
236 std::next(coordinate_dofs.begin(), i * gdim));
239 auto coeff_cell = coeffs.row(cell);
240 fn(&value, coeff_cell.data(), constant_values.data(),
241 coordinate_dofs.data(), &local_facet,
242 &perms[cell * facets.size() + local_facet], cell_info[cell]);
248 template <
typename T>
249 T assemble_interior_facets(
250 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
251 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
252 const std::uint8_t*,
const std::uint32_t)>& fn,
253 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
255 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
256 const std::vector<std::uint32_t>& cell_info,
257 const std::vector<std::uint8_t>& perms)
259 const int gdim = mesh.geometry().dim();
260 const int tdim = mesh.topology().dim();
263 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
266 const int num_dofs_g = x_dofmap.num_links(0);
267 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
268 = mesh.geometry().x();
271 std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
272 std::vector<T> coeff_array(2 * offsets.back());
273 assert(offsets.back() == coeffs.cols());
275 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
277 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
282 const int offset_g = gdim * num_dofs_g;
283 for (std::int32_t f : active_facets)
286 auto cells = f_to_c->links(f);
287 assert(
cells.size() == 2);
289 const int facets_per_cell = c_to_f->num_links(cells[0]);
292 std::array<int, 2> local_facet;
293 for (
int i = 0; i < 2; ++i)
295 auto facets = c_to_f->links(cells[i]);
296 auto it = std::find(facets.begin(), facets.end(), f);
297 assert(it != facets.end());
298 local_facet[i] = std::distance(facets.begin(), it);
302 auto x_dofs0 = x_dofmap.links(cells[0]);
303 auto x_dofs1 = x_dofmap.links(cells[1]);
304 for (
int i = 0; i < num_dofs_g; ++i)
306 for (
int j = 0; j < gdim; ++j)
308 coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
309 coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
315 auto coeff_cell0 = coeffs.row(cells[0]);
316 auto coeff_cell1 = coeffs.row(cells[1]);
319 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
322 const int num_entries = offsets[i + 1] - offsets[i];
323 std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
324 std::next(coeff_array.begin(), 2 * offsets[i]));
325 std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
326 std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
329 const std::array perm{perms[
cells[0] * facets_per_cell + local_facet[0]],
330 perms[
cells[1] * facets_per_cell + local_facet[1]]};
331 fn(&value, coeff_array.data(), constant_values.data(),
332 coordinate_dofs.data(), local_facet.data(), perm.data(),
333 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
Eigen::Array< typename U::scalar_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:400