11 #include <Eigen/Dense>
12 #include <dolfinx/common/IndexMap.h>
13 #include <dolfinx/common/types.h>
14 #include <dolfinx/function/Constant.h>
15 #include <dolfinx/function/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 Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& 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 Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
51 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& 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 Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
63 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& 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 function::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.data(),
85 array.data() + array.size());
89 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
92 const bool needs_permutation_data = M.needs_permutation_data();
93 if (needs_permutation_data)
94 mesh->topology_mutable().create_entity_permutations();
95 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
96 = needs_permutation_data
97 ? mesh->topology().get_cell_permutation_info()
98 : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
101 for (
int i : M.integral_ids(IntegralType::cell))
103 const auto& fn = M.kernel(IntegralType::cell, i);
104 const std::vector<std::int32_t>& active_cells
105 = M.domains(IntegralType::cell, i);
106 value += fem::impl::assemble_cells(mesh->geometry(), active_cells, fn,
107 coeffs, constant_values, cell_info);
110 if (M.num_integrals(IntegralType::exterior_facet) > 0
111 or M.num_integrals(IntegralType::interior_facet) > 0)
114 mesh->topology_mutable().create_entities(tdim - 1);
115 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
117 const int facets_per_cell
119 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
120 = needs_permutation_data
121 ? mesh->topology().get_facet_permutations()
122 : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
123 facets_per_cell, num_cells);
125 for (
int i : M.integral_ids(IntegralType::exterior_facet))
127 const auto& fn = M.kernel(IntegralType::exterior_facet, i);
128 const std::vector<std::int32_t>& active_facets
129 = M.domains(IntegralType::exterior_facet, i);
130 value += fem::impl::assemble_exterior_facets(
131 *mesh, active_facets, fn, coeffs, constant_values, cell_info, perms);
134 const std::vector<int> c_offsets = M.coefficient_offsets();
135 for (
int i : M.integral_ids(IntegralType::interior_facet))
137 const auto& fn = M.kernel(IntegralType::interior_facet, i);
138 const std::vector<std::int32_t>& active_facets
139 = M.domains(IntegralType::interior_facet, i);
140 value += fem::impl::assemble_interior_facets(
141 *mesh, active_facets, fn, coeffs, c_offsets, constant_values,
149 template <
typename T>
151 const mesh::Geometry& geometry,
152 const std::vector<std::int32_t>& active_cells,
153 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
154 const std::uint8_t*,
const std::uint32_t)>& fn,
155 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
157 const std::vector<T>& constant_values,
158 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info)
160 const int gdim = geometry.dim();
163 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
166 const int num_dofs_g = x_dofmap.num_links(0);
167 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
171 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
172 coordinate_dofs(num_dofs_g, gdim);
176 for (std::int32_t c : active_cells)
179 auto x_dofs = x_dofmap.links(c);
180 for (
int i = 0; i < num_dofs_g; ++i)
181 for (
int j = 0; j < gdim; ++j)
182 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
184 auto coeff_cell = coeffs.row(c);
185 fn(&value, coeff_cell.data(), constant_values.data(),
186 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
192 template <
typename T>
193 T assemble_exterior_facets(
194 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
195 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
196 const std::uint8_t*,
const std::uint32_t)>& fn,
197 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
199 const std::vector<T>& constant_values,
200 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
201 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
203 const int gdim = mesh.geometry().dim();
204 const int tdim = mesh.topology().dim();
207 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
210 const int num_dofs_g = x_dofmap.num_links(0);
211 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
212 = mesh.geometry().x();
215 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
216 coordinate_dofs(num_dofs_g, gdim);
218 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
220 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
225 for (std::int32_t facet : active_facets)
228 assert(f_to_c->num_links(facet) == 1);
229 const int cell = f_to_c->links(facet)[0];
232 auto facets = c_to_f->links(cell);
234 = std::find(facets.data(), facets.data() + facets.rows(), facet);
235 assert(it != (facets.data() + facets.rows()));
236 const int local_facet = std::distance(facets.data(), it);
238 auto x_dofs = x_dofmap.links(cell);
239 for (
int i = 0; i < num_dofs_g; ++i)
240 for (
int j = 0; j < gdim; ++j)
241 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
243 auto coeff_cell = coeffs.row(cell);
244 fn(&value, coeff_cell.data(), constant_values.data(),
245 coordinate_dofs.data(), &local_facet, &perms(local_facet, cell),
252 template <
typename T>
253 T assemble_interior_facets(
254 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
255 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
256 const std::uint8_t*,
const std::uint32_t)>& fn,
257 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
259 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
260 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
261 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
263 const int gdim = mesh.geometry().dim();
264 const int tdim = mesh.topology().dim();
267 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
270 const int num_dofs_g = x_dofmap.num_links(0);
271 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
272 = mesh.geometry().x();
275 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
276 coordinate_dofs(2 * num_dofs_g, gdim);
277 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
278 assert(offsets.back() == coeffs.cols());
280 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
282 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
287 for (std::int32_t f : active_facets)
290 auto cells = f_to_c->links(f);
291 assert(
cells.rows() == 2);
294 std::array<int, 2> local_facet;
295 for (
int i = 0; i < 2; ++i)
297 auto facets = c_to_f->links(cells[i]);
299 = std::find(facets.data(), facets.data() + facets.rows(), f);
300 assert(it != (facets.data() + facets.rows()));
301 local_facet[i] = std::distance(facets.data(), it);
305 auto x_dofs0 = x_dofmap.links(cells[0]);
306 auto x_dofs1 = x_dofmap.links(cells[1]);
307 for (
int i = 0; i < num_dofs_g; ++i)
309 for (
int j = 0; j < gdim; ++j)
311 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
312 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
318 auto coeff_cell0 = coeffs.row(cells[0]);
319 auto coeff_cell1 = coeffs.row(cells[1]);
322 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
325 const int num_entries = offsets[i + 1] - offsets[i];
326 coeff_array.segment(2 * offsets[i], num_entries)
327 = coeff_cell0.segment(offsets[i], num_entries);
328 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
329 = coeff_cell1.segment(offsets[i], num_entries);
332 const std::array perm{perms(local_facet[0], cells[0]),
333 perms(local_facet[1], cells[1])};
334 fn(&value, coeff_array.data(), constant_values.data(),
335 coordinate_dofs.data(), local_facet.data(), perm.data(),
336 cell_info[cells[0]]);
void cells(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const fem::DofMap *, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: SparsityPatternBuilder.cpp:18
T assemble_scalar(const Form< T > &M)
Assemble functional into scalar. Caller is responsible for accumulation across processes.
Definition: assembler.h:38
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:407
int cell_num_entities(mesh::CellType type, int dim)
Number of entities of dimension dim.
Definition: cell_types.cpp:381