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,
43 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
46 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
47 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
48 const std::uint8_t*,
const std::uint32_t)>& kernel,
49 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
51 const Eigen::Array<T, Eigen::Dynamic, 1>& constants);
55 void assemble_exterior_facets(
56 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
57 const std::int32_t*,
const T*)>& mat_set_values,
58 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
61 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
62 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
63 const std::uint8_t*,
const std::uint32_t)>& fn,
64 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
66 const Eigen::Array<T, Eigen::Dynamic, 1> constants);
70 void assemble_interior_facets(
71 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
72 const std::int32_t*,
const T*)>& mat_set_values,
73 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
74 const DofMap& dofmap0,
const DofMap& dofmap1,
const std::vector<bool>& bc0,
75 const std::vector<bool>& bc1,
76 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
77 const std::uint8_t*,
const std::uint32_t)>& kernel,
78 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
80 const std::vector<int>& offsets,
81 const Eigen::Array<T, Eigen::Dynamic, 1>& constants);
86 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
87 const std::int32_t*,
const T*)>& mat_set_values,
88 const Form<T>& a,
const std::vector<bool>& bc0,
89 const std::vector<bool>& bc1)
91 std::shared_ptr<const mesh::Mesh> mesh = a.
mesh();
95 std::shared_ptr<const fem::DofMap> dofmap0 = a.
function_space(0)->dofmap();
96 std::shared_ptr<const fem::DofMap> dofmap1 = a.
function_space(1)->dofmap();
104 throw std::runtime_error(
"Unset constant in Form");
105 const Eigen::Array<T, Eigen::Dynamic, 1> constants =
pack_constants(a);
108 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
112 for (
int i = 0; i < integrals.
num_integrals(IntegralType::cell); ++i)
115 const std::vector<std::int32_t>& active_cells
117 fem::impl::assemble_cells<T>(mat_set_values, *mesh, active_cells, dofs0,
118 dofs1, bc0, bc1, fn, coeffs, constants);
121 for (
int i = 0; i < integrals.
num_integrals(IntegralType::exterior_facet);
126 const std::vector<std::int32_t>& active_facets
128 fem::impl::assemble_exterior_facets<T>(mat_set_values, *mesh, active_facets,
129 dofs0, dofs1, bc0, bc1, fn, coeffs,
133 const std::vector<int> c_offsets = a.
coefficients().offsets();
134 for (
int i = 0; i < integrals.
num_integrals(IntegralType::interior_facet);
139 const std::vector<std::int32_t>& active_facets
141 fem::impl::assemble_interior_facets<T>(mat_set_values, *mesh, active_facets,
142 *dofmap0, *dofmap1, bc0, bc1, fn,
143 coeffs, c_offsets, constants);
147 template <
typename T>
149 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
150 const std::int32_t*,
const T*)>& mat_set,
151 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
154 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
155 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
156 const std::uint8_t*,
const std::uint32_t)>& kernel,
157 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
159 const Eigen::Array<T, Eigen::Dynamic, 1>& constants)
168 const int num_dofs_g = x_dofmap.
num_links(0);
169 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
173 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
174 coordinate_dofs(num_dofs_g, gdim);
175 const int num_dofs0 = dofmap0.
links(0).size();
176 const int num_dofs1 = dofmap1.
links(0).size();
177 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae(
178 num_dofs0, num_dofs1);
180 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
184 for (std::int32_t c : active_cells)
187 auto x_dofs = x_dofmap.
links(c);
188 for (
int i = 0; i < x_dofs.rows(); ++i)
189 for (
int j = 0; j < gdim; ++j)
190 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
193 std::fill(Ae.data(), Ae.data() + num_dofs0 * num_dofs1, 0);
194 kernel(Ae.data(), coeffs.row(c).data(), constants.data(),
195 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
198 auto dofs0 = dofmap0.
links(c);
199 auto dofs1 = dofmap1.
links(c);
202 for (Eigen::Index i = 0; i < Ae.rows(); ++i)
210 for (Eigen::Index j = 0; j < Ae.cols(); ++j)
217 mat_set(dofs0.size(), dofs0.data(), dofs1.size(), dofs1.data(), Ae.data());
221 template <
typename T>
222 void assemble_exterior_facets(
223 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
224 const std::int32_t*,
const T*)>& mat_set_values,
225 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
228 const std::vector<bool>& bc0,
const std::vector<bool>& bc1,
229 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
230 const std::uint8_t*,
const std::uint32_t)>& kernel,
231 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
233 const Eigen::Array<T, Eigen::Dynamic, 1> constants)
247 const int num_dofs_g = x_dofmap.
num_links(0);
248 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
252 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
253 coordinate_dofs(num_dofs_g, gdim);
254 const int num_dofs0 = dofmap0.
links(0).size();
255 const int num_dofs1 = dofmap1.
links(0).size();
256 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae(
257 num_dofs0, num_dofs1);
259 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
261 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
269 for (std::int32_t f : active_facets)
271 auto cells = f_to_c->links(f);
272 assert(cells.rows() == 1);
275 auto facets = c_to_f->links(cells[0]);
276 const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
277 assert(it != (facets.data() + facets.rows()));
278 const int local_facet = std::distance(facets.data(), it);
281 auto x_dofs = x_dofmap.
links(cells[0]);
282 for (
int i = 0; i < num_dofs_g; ++i)
283 for (
int j = 0; j < gdim; ++j)
284 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
287 const std::uint8_t perm = perms(local_facet, cells[0]);
288 std::fill(Ae.data(), Ae.data() + num_dofs0 * num_dofs1, 0);
289 kernel(Ae.data(), coeffs.row(cells[0]).data(), constants.data(),
290 coordinate_dofs.data(), &local_facet, &perm, cell_info[cells[0]]);
293 auto dmap0 = dofmap0.
links(cells[0]);
294 auto dmap1 = dofmap1.
links(cells[0]);
297 for (Eigen::Index i = 0; i < Ae.rows(); ++i)
305 for (Eigen::Index j = 0; j < Ae.cols(); ++j)
312 mat_set_values(dmap0.size(), dmap0.data(), dmap1.size(), dmap1.data(),
317 template <
typename T>
318 void assemble_interior_facets(
319 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
320 const std::int32_t*,
const T*)>& mat_set_values,
321 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
322 const DofMap& dofmap0,
const DofMap& dofmap1,
const std::vector<bool>& bc0,
323 const std::vector<bool>& bc1,
324 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
325 const std::uint8_t*,
const std::uint32_t)>& fn,
326 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
328 const std::vector<int>& offsets,
329 const Eigen::Array<T, Eigen::Dynamic, 1>& constants)
343 const int num_dofs_g = x_dofmap.
num_links(0);
345 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
349 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
350 coordinate_dofs(2 * num_dofs_g, gdim);
351 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
352 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
353 assert(offsets.back() == coeffs.cols());
356 Eigen::Array<std::int32_t, Eigen::Dynamic, 1> dmapjoint0, dmapjoint1;
358 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
360 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
368 for (std::int32_t facet_index : active_facets)
371 auto cells = c->links(facet_index);
372 assert(cells.rows() == 2);
375 auto facets0 = c_to_f->links(cells[0]);
376 const auto* it0 = std::find(facets0.data(), facets0.data() + facets0.rows(),
378 assert(it0 != (facets0.data() + facets0.rows()));
379 const int local_facet0 = std::distance(facets0.data(), it0);
380 auto facets1 = c_to_f->links(cells[1]);
381 const auto* it1 = std::find(facets1.data(), facets1.data() + facets1.rows(),
383 assert(it1 != (facets1.data() + facets1.rows()));
384 const int local_facet1 = std::distance(facets1.data(), it1);
386 const std::array local_facet{local_facet0, local_facet1};
387 const std::array perm{perms(local_facet[0], cells[0]),
388 perms(local_facet[1], cells[1])};
391 auto x_dofs0 = x_dofmap.
links(cells[0]);
392 auto x_dofs1 = x_dofmap.
links(cells[1]);
393 for (
int i = 0; i < num_dofs_g; ++i)
395 for (
int j = 0; j < gdim; ++j)
397 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
398 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
403 auto dmap0_cell0 = dofmap0.
cell_dofs(cells[0]);
404 auto dmap0_cell1 = dofmap0.
cell_dofs(cells[1]);
405 dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
406 dmapjoint0.head(dmap0_cell0.size()) = dmap0_cell0;
407 dmapjoint0.tail(dmap0_cell1.size()) = dmap0_cell1;
409 auto dmap1_cell0 = dofmap1.
cell_dofs(cells[0]);
410 auto dmap1_cell1 = dofmap1.
cell_dofs(cells[1]);
411 dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
412 dmapjoint1.head(dmap1_cell0.size()) = dmap1_cell0;
413 dmapjoint1.tail(dmap1_cell1.size()) = dmap1_cell1;
417 auto coeff_cell0 = coeffs.row(cells[0]);
418 auto coeff_cell1 = coeffs.row(cells[1]);
421 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
424 const int num_entries = offsets[i + 1] - offsets[i];
425 coeff_array.segment(2 * offsets[i], num_entries)
426 = coeff_cell0.segment(offsets[i], num_entries);
427 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
428 = coeff_cell1.segment(offsets[i], num_entries);
432 Ae.setZero(dmapjoint0.size(), dmapjoint1.size());
433 fn(Ae.data(), coeff_array.data(), constants.data(), coordinate_dofs.data(),
434 local_facet.data(), perm.data(), cell_info[cells[0]]);
439 for (Eigen::Index i = 0; i < dmapjoint0.size(); ++i)
441 if (bc0[dmapjoint0[i]])
447 for (Eigen::Index j = 0; j < dmapjoint1.size(); ++j)
449 if (bc1[dmapjoint1[j]])
454 mat_set_values(dmapjoint0.size(), dmapjoint0.data(), dmapjoint1.size(),
455 dmapjoint1.data(), Ae.data());