12 #include <Eigen/Dense>
13 #include <dolfinx/function/FunctionSpace.h>
14 #include <dolfinx/graph/AdjacencyList.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>
22 namespace dolfinx::fem::impl
33 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
34 const std::int32_t*,
const T*)>& mat_set_values,
35 const Form<T>& a,
const std::vector<bool>& bc0,
36 const std::vector<bool>& bc1);
41 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
42 const std::int32_t*,
const T*)>& mat_set_values,
44 const std::vector<std::int32_t>& active_cells,
47 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
48 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
49 const std::uint8_t*,
const std::uint32_t)>& kernel,
50 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
52 const Eigen::Array<T, Eigen::Dynamic, 1>& constants,
53 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info);
57 void assemble_exterior_facets(
58 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
59 const std::int32_t*,
const T*)>& mat_set_values,
60 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
63 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
64 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
65 const std::uint8_t*,
const std::uint32_t)>& fn,
66 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
68 const Eigen::Array<T, Eigen::Dynamic, 1> constants,
69 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
70 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms);
74 void assemble_interior_facets(
75 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
76 const std::int32_t*,
const T*)>& mat_set_values,
77 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
78 const DofMap& dofmap0,
const DofMap& dofmap1,
const std::vector<bool>& bc0,
79 const std::vector<bool>& bc1,
80 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
81 const std::uint8_t*,
const std::uint32_t)>& kernel,
82 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
84 const std::vector<int>& offsets,
85 const Eigen::Array<T, Eigen::Dynamic, 1>& constants,
86 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
87 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms);
92 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
93 const std::int32_t*,
const T*)>& mat_set_values,
94 const Form<T>& a,
const std::vector<bool>& bc0,
95 const std::vector<bool>& bc1)
97 std::shared_ptr<const mesh::Mesh> mesh = a.
mesh();
99 const int tdim = mesh->topology().dim();
100 const std::int32_t num_cells
101 = mesh->topology().connectivity(tdim, 0)->num_nodes();
104 std::shared_ptr<const fem::DofMap> dofmap0
106 std::shared_ptr<const fem::DofMap> dofmap1
114 const Eigen::Array<T, Eigen::Dynamic, 1> constants =
pack_constants(a);
117 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
121 if (needs_permutation_data)
122 mesh->topology_mutable().create_entity_permutations();
123 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
124 = needs_permutation_data
125 ? mesh->topology().get_cell_permutation_info()
126 : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
130 const auto& fn = a.
kernel(IntegralType::cell, i);
131 const std::vector<std::int32_t>& active_cells
132 = a.
domains(IntegralType::cell, i);
133 impl::assemble_cells<T>(mat_set_values, mesh->geometry(), active_cells,
134 dofs0, dofs1, bc0, bc1, fn, coeffs, constants,
142 mesh->topology_mutable().create_entities(tdim - 1);
143 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
145 const int facets_per_cell
147 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
148 = needs_permutation_data
149 ? mesh->topology().get_facet_permutations()
150 : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
151 facets_per_cell, num_cells);
153 for (
int i : a.
integral_ids(IntegralType::exterior_facet))
155 const auto& fn = a.
kernel(IntegralType::exterior_facet, i);
156 const std::vector<std::int32_t>& active_facets
157 = a.
domains(IntegralType::exterior_facet, i);
158 impl::assemble_exterior_facets<T>(mat_set_values, *mesh, active_facets,
159 dofs0, dofs1, bc0, bc1, fn, coeffs,
160 constants, cell_info, perms);
164 for (
int i : a.
integral_ids(IntegralType::interior_facet))
166 const auto& fn = a.
kernel(IntegralType::interior_facet, i);
167 const std::vector<std::int32_t>& active_facets
168 = a.
domains(IntegralType::interior_facet, i);
169 impl::assemble_interior_facets<T>(
170 mat_set_values, *mesh, active_facets, *dofmap0, *dofmap1, bc0, bc1,
171 fn, coeffs, c_offsets, constants, cell_info, perms);
176 template <
typename T>
178 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
179 const std::int32_t*,
const T*)>& mat_set,
181 const std::vector<std::int32_t>& active_cells,
184 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
185 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
186 const std::uint8_t*,
const std::uint32_t)>& kernel,
187 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
189 const Eigen::Array<T, Eigen::Dynamic, 1>& constants,
190 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info)
192 const int gdim = geometry.
dim();
198 const int num_dofs_g = x_dofmap.
num_links(0);
199 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
203 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
204 coordinate_dofs(num_dofs_g, gdim);
205 const int num_dofs0 = dofmap0.
links(0).size();
206 const int num_dofs1 = dofmap1.
links(0).size();
207 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae(
208 num_dofs0, num_dofs1);
211 for (std::int32_t c : active_cells)
214 auto x_dofs = x_dofmap.
links(c);
215 for (
int i = 0; i < x_dofs.rows(); ++i)
216 for (
int j = 0; j < gdim; ++j)
217 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
220 std::fill(Ae.data(), Ae.data() + num_dofs0 * num_dofs1, 0);
221 kernel(Ae.data(), coeffs.row(c).data(), constants.data(),
222 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
225 auto dofs0 = dofmap0.
links(c);
226 auto dofs1 = dofmap1.
links(c);
229 for (Eigen::Index i = 0; i < Ae.rows(); ++i)
237 for (Eigen::Index j = 0; j < Ae.cols(); ++j)
244 mat_set(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(), Ae.data());
248 template <
typename T>
249 void assemble_exterior_facets(
250 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
251 const std::int32_t*,
const T*)>& mat_set_values,
252 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
255 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
256 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
257 const std::uint8_t*,
const std::uint32_t)>& kernel,
258 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
260 const Eigen::Array<T, Eigen::Dynamic, 1> constants,
261 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
262 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
271 const int num_dofs_g = x_dofmap.
num_links(0);
272 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
276 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
277 coordinate_dofs(num_dofs_g, gdim);
278 const int num_dofs0 = dofmap0.
links(0).size();
279 const int num_dofs1 = dofmap1.
links(0).size();
280 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae(
281 num_dofs0, num_dofs1);
288 for (std::int32_t f : active_facets)
290 auto cells = f_to_c->links(f);
291 assert(cells.rows() == 1);
294 auto facets = c_to_f->links(cells[0]);
295 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
296 assert(it != (facets.data() + facets.rows()));
297 const int local_facet = std::distance(facets.data(), it);
300 auto x_dofs = x_dofmap.
links(cells[0]);
301 for (
int i = 0; i < num_dofs_g; ++i)
302 for (
int j = 0; j < gdim; ++j)
303 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
306 std::fill(Ae.data(), Ae.data() + num_dofs0 * num_dofs1, 0);
307 kernel(Ae.data(), coeffs.row(cells[0]).data(), constants.data(),
308 coordinate_dofs.data(), &local_facet, &perms(local_facet, cells[0]),
309 cell_info[cells[0]]);
312 auto dmap0 = dofmap0.
links(cells[0]);
313 auto dmap1 = dofmap1.
links(cells[0]);
316 for (Eigen::Index i = 0; i < Ae.rows(); ++i)
324 for (Eigen::Index j = 0; j < Ae.cols(); ++j)
331 mat_set_values(dmap0.size(), dmap0.data(), dmap1.size(), dmap1.data(),
336 template <
typename T>
337 void assemble_interior_facets(
338 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
339 const std::int32_t*,
const T*)>& mat_set_values,
340 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
341 const DofMap& dofmap0,
const DofMap& dofmap1,
const std::vector<bool>& bc0,
342 const std::vector<bool>& bc1,
343 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
344 const std::uint8_t*,
const std::uint32_t)>& fn,
345 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
347 const std::vector<int>& offsets,
348 const Eigen::Array<T, Eigen::Dynamic, 1>& constants,
349 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
350 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
359 const int num_dofs_g = x_dofmap.
num_links(0);
361 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
365 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
366 coordinate_dofs(2 * num_dofs_g, gdim);
367 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
368 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
369 assert(offsets.back() == coeffs.cols());
372 Eigen::Array<std::int32_t, Eigen::Dynamic, 1> dmapjoint0, dmapjoint1;
379 for (std::int32_t facet_index : active_facets)
382 auto cells = c->links(facet_index);
383 assert(cells.rows() == 2);
386 auto facets0 = c_to_f->links(cells[0]);
387 const auto* it0 = std::find(facets0.data(), facets0.data() + facets0.rows(),
389 assert(it0 != (facets0.data() + facets0.rows()));
390 const int local_facet0 = std::distance(facets0.data(), it0);
391 auto facets1 = c_to_f->links(cells[1]);
392 const auto* it1 = std::find(facets1.data(), facets1.data() + facets1.rows(),
394 assert(it1 != (facets1.data() + facets1.rows()));
395 const int local_facet1 = std::distance(facets1.data(), it1);
397 const std::array local_facet{local_facet0, local_facet1};
400 auto x_dofs0 = x_dofmap.
links(cells[0]);
401 auto x_dofs1 = x_dofmap.
links(cells[1]);
402 for (
int i = 0; i < num_dofs_g; ++i)
404 for (
int j = 0; j < gdim; ++j)
406 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
407 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
412 auto dmap0_cell0 = dofmap0.
cell_dofs(cells[0]);
413 auto dmap0_cell1 = dofmap0.
cell_dofs(cells[1]);
414 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
415 dmapjoint0.head(dmap0_cell0.size()) = dmap0_cell0;
416 dmapjoint0.tail(dmap0_cell1.size()) = dmap0_cell1;
418 auto dmap1_cell0 = dofmap1.
cell_dofs(cells[0]);
419 auto dmap1_cell1 = dofmap1.
cell_dofs(cells[1]);
420 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
421 dmapjoint1.head(dmap1_cell0.size()) = dmap1_cell0;
422 dmapjoint1.tail(dmap1_cell1.size()) = dmap1_cell1;
426 auto coeff_cell0 = coeffs.row(cells[0]);
427 auto coeff_cell1 = coeffs.row(cells[1]);
430 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
433 const int num_entries = offsets[i + 1] - offsets[i];
434 coeff_array.segment(2 * offsets[i], num_entries)
435 = coeff_cell0.segment(offsets[i], num_entries);
436 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
437 = coeff_cell1.segment(offsets[i], num_entries);
441 Ae.setZero(dmapjoint0.size(), dmapjoint1.size());
442 const std::array perm{perms(local_facet[0], cells[0]),
443 perms(local_facet[1], cells[1])};
444 fn(Ae.data(), coeff_array.data(), constants.data(), coordinate_dofs.data(),
445 local_facet.data(), perm.data(), cell_info[cells[0]]);
450 for (Eigen::Index i = 0; i < dmapjoint0.size(); ++i)
452 if (bc0[dmapjoint0[i]])
458 for (Eigen::Index j = 0; j < dmapjoint1.size(); ++j)
460 if (bc1[dmapjoint1[j]])
465 mat_set_values(dmapjoint0.size(), dmapjoint0.data(), dmapjoint1.size(),
466 dmapjoint1.data(), Ae.data());